#include "Cluster.hpp"
#include "GraphHandler.hpp"
#include "GraphPredicates.hpp"
#include "IsoGraphEditDistance.hpp"
#include "Partition.h"
#include "SetPartition.hpp"
#include "GraphExploration.hpp"
#include "EditOperationMap.hpp"

#include <sys/types.h>
#include <dirent.h>
#include <errno.h>
#include <sys/types.h>
#include <sys/stat.h>
#include <unistd.h>
#include <sstream>
#include <fstream>
#include <iostream>
#include <boost/algorithm/string.hpp>
#include <boost/math/special_functions/binomial.hpp>
#include <boost/graph/adj_list_serialize.hpp>
#include <boost/graph/connected_components.hpp>
#include <boost/timer.hpp>

// include headers that implement an archive
#include <boost/archive/binary_oarchive.hpp>
#include <boost/archive/binary_iarchive.hpp>




// To serialize std::list
#include <boost/serialization/list.hpp>
#include <boost/serialization/vector.hpp>
#include <boost/serialization/map.hpp>
#include <boost/serialization/set.hpp>


//----------------------------------------------------------------
GraphHandler::GraphHandler()
{
    m_bMetaFiltered = false;

    m_iIDCounter = 0;

}
//----------------------------------------------------------------
GraphHandler::~GraphHandler()
{
}
//----------------------------------------------------------------
// Returns (name, type)
pair<string,string> GraphHandler::GetNameTypeAttr
        (const property_tree::ptree::value_type v) const
{
    using boost::property_tree::ptree;
    string sName, sType;

    BOOST_FOREACH(const ptree::value_type &i,
    v.second.get_child("<xmlattr>"))
    {
        if( ((string)i.first.data()).compare("name") == 0 )
            sName = i.second.data();
        else
            sType = i.second.data();
    }

    return make_pair(sName, sType);
}
//----------------------------------------------------------------
//template <typename GeneralGraph>
//pair< bool, typename graph_traits<GeneralGraph>::vertex_iterator >
pair< bool, GraphHandler::Iter_t >
    GraphHandler::FindVertex(string sName, const Graph& g) const
{

    //typedef pair<bool, typename graph_traits<GeneralGraph>::vertex_iterator> pair_t;
    //typedef typename graph_traits<GeneralGraph>::vertex_iterator vertex_iter;
    typedef pair< bool, Iter_t > pair_t;

    //std::pair<vertex_iter, vertex_iter> vp;
    std::pair<Iter_t, Iter_t> vp;
    for (vp = vertices(g); vp.first != vp.second; ++vp.first)
    {
        //cout << "Comparing:" << g[*vp.first].sName << endl;
        //if( g[*vp.first].sName.compare(sName) == 0 )
        //const string stmp = m_IDTypes.find(get(vertex_name, g, *vp.first));
        map<int, string>::const_iterator it = m_IDTypes.find(get(vertex_name, g, *vp.first));
        //string stmp = (*it);
        if( it->second.compare(sName) == 0 )
            return pair_t(true, vp.first);
    }
    return pair_t(false, vp.first);
}
//----------------------------------------------------------------
bool GraphHandler::LoadSingleGraph(string sFilename)
{
    // Create an empty property tree object
    using boost::property_tree::ptree;
    ptree pt;

    // This maps from name -> type
    typedef map< string, string > map_t;
    map_t nameData;

    // This is used for collecting meta-data, we should only count
    // each category once per graph. This is used to check that.
    map<string, bool> TypeCountedMap;



    // A temporary storage for the Graph being created
    Graph kTmpG;


    // A local map for this specific graph, from type->vertex (in local graph kTmpG)
    map<string, Vertex> TypeVertexMap;

    // Load the XML file into the property tree. If reading fails
    // (cannot open file, parse error), an exception is thrown.
    read_xml(sFilename, pt);
    //get(graph_name, kTmpG) = sFilename;

    size_t found;
    found=sFilename.find_last_of("//");

    //sFilename.find_last_of("//")
    set_property(kTmpG, boost::graph_name, sFilename.substr(found+1));

    // We shall first map the different names to types, create a map from Name -> Type
    // since we're not interested in specific rooms, we want the general topology.

    // Go through and pick out the "space" categories from the file
    BOOST_FOREACH(ptree::value_type &v,
            pt.get_child("MITquest"))
    {
            if( ((string)v.first.data()).compare("space") == 0 )
            {
                // This is a "room", ie a semantic category.
                // Get the xmlattribute type.
                pair<string,string> pTmp = GetNameTypeAttr(v);


               pair<set<string>::iterator,bool> kRet = m_Types.insert(pTmp.second);
               if( kRet.second )
               {
                    // Indeed a new vertex found, add it.
                    MetaVertex pNew;
                    pNew = add_vertex(m_MetaGraph);
                    m_MetaGraph[pNew].sName = pTmp.second;
                    m_MetaGraph[pNew].iFreq = 0;
                    // Record which vertex this name belongs to.
                    m_TypeVertexMap[pTmp.second] = pNew;

                    m_TypesID[pTmp.second] = m_iIDCounter;
                    m_IDTypes[m_iIDCounter] = pTmp.second;
                    m_iIDCounter++;

                    put(vertex_name, m_MetaGraph, pNew, m_TypesID[pTmp.second]);
                    //TypeCountedMap[pTmp.second] = true;
               }

                if( TypeCountedMap.find(pTmp.second) == TypeCountedMap.end() )
                {
                    // We should only increment the meta freq counter once per graph
                    m_MetaGraph[m_TypeVertexMap[pTmp.second]].iFreq++;

                    // Add the new vertex to the local graph.
                    Vertex pNew;
                    pNew = add_vertex(kTmpG);
                    //kTmpG[pNew].sName = pTmp.second;
                    //put(property, g, x, v) 	Set property value for vertex or edge x to v.
                    put(vertex_name, kTmpG, pNew, m_TypesID[pTmp.second]);
                    //put(vertex_id, kTmpG, pNew, m_TypesID[pTmp.second]);
                    //cout << "Adding: " << kTmpG[pNew].sName << endl;
                    TypeVertexMap[pTmp.second] = pNew;

                    TypeCountedMap[pTmp.second] = true;
                }
               //kTmpG

               // Associate this name with the category.
               nameData[pTmp.first] = pTmp.second;
            }
        }


    // Now construct the graph, go through all the rooms again
    BOOST_FOREACH(ptree::value_type &v,
            pt.get_child("MITquest"))
    {
            if( ((string)v.first.data()).compare("space") == 0 )
            {
                pair<string,string> pOrigin = GetNameTypeAttr(v);
                string sOrigin = pOrigin.second;
                // pOrigin.second has the vertex from which we are connecting

                // Get the vertex which corresponds to this name
                MetaVertex vSource = m_TypeVertexMap[sOrigin];


                BOOST_FOREACH(ptree::value_type &i,
                        v.second)
                {
                    // Compare if its portal

                    if( ((string)i.first.data()).compare("portal") == 0 )
                    {


                        BOOST_FOREACH(ptree::value_type &j,
                                i.second.get_child("<xmlattr>"))
                                {
                                    if( ((string)j.first.data()).compare("target") == 0 )
                                    {
                                        //string sTarget = ((string)j.second.data());
                                        // Now ensure that this specific room is in our mapping
                                        // from specific rooms -> types. (this is due to the data format,
                                        // some of the sources does not seem to be valid rooms)
                                        // If so, add an edge to the topology and update meta-topology.
                                        map<string,string>::const_iterator itTmp =
                                            nameData.find(((string)j.second.data()));



                                        if( itTmp != nameData.end()
                                            && m_Types.find(itTmp->second) != m_Types.end() )
                                        {
                                                string sTarget = itTmp->second;

                                                pair<Edge, bool> kRet2 = edge(TypeVertexMap[sOrigin], TypeVertexMap[sTarget], kTmpG);

                                                // Don't allow self-edges
                                                if( !kRet2.second && TypeVertexMap[sOrigin] != TypeVertexMap[sTarget] )
                                                // Add edge in local graph
                                                add_edge(TypeVertexMap[sOrigin],TypeVertexMap[sTarget], kTmpG);
                                                // This is a valid connection, add it.
                                                //cout << "Target: " << sTarget << endl;
                                                MetaVertex vTarget = m_TypeVertexMap[sTarget];

                                                // Add the edge if it does not already exist.
                                                // edge(u,v,g) returns pair<edge_descriptor, bool>
                                                pair<MetaEdge, bool> kRet = edge( vSource, vTarget, m_MetaGraph);

                                                if( !kRet.second )
                                                {
                                                    // Add the edge, it does not exist
                                                    add_edge(vSource, vTarget,
                                                    EdgeMetaProperties(1), m_MetaGraph);
                                                    //m_MetaGraph
                                                }
                                                else m_MetaGraph[kRet.first].iFreq++;
                                        }
                                    }
                                }

                    }



                }
            }
    }

    m_Graphs.push_back(kTmpG);

    return true;

}
//----------------------------------------------------------------
int GraphHandler::CountFilesInDir(string sDir) const
{
    DIR *dp;
    stringstream fLoad;
    struct dirent *dirp;

    if((dp  = opendir(sDir.c_str())) == NULL)
    {   cout << "Error(" << errno << ") opening " << sDir << endl;
        return 0;
    }
    int iX = 0;
    while ((dirp = readdir(dp)) != NULL ){
        if( strlen(dirp->d_name) != 0 && dirp->d_name[0] != '.' )
        {
            fLoad << sDir << "/" << dirp->d_name;
            struct stat statTmp;
            stat(fLoad.str().c_str(), &statTmp);

            // Check if it is a file or directory
            if( !S_ISDIR(statTmp.st_mode) )
                iX++;
            else
            {   iX += CountFilesInDir(string(fLoad.str()));
                fLoad.str(""); }
        }
    }
    return iX;
}
//----------------------------------------------------------------
int GraphHandler::LoadAllGraphs(string sDir, int iLimit)
{
    DIR *dp;
    stringstream fLoad;
    struct dirent *dirp;
    if((dp  = opendir(sDir.c_str())) == NULL)
    {
        cout << "Error(" << errno << ") opening " << sDir << endl;
        return 0;
    }

    int iCounter = 0;

    int iFiles = CountFilesInDir(sDir);
    m_Graphs.reserve(iFiles);

    // Load all the files (except current directory(.) and previous dir(..))
    while ((dirp = readdir(dp)) != NULL )
    {
        if( strlen(dirp->d_name) != 0 && dirp->d_name[0] != '.' )
        {
            fLoad << sDir << "/" << dirp->d_name;

            struct stat statTmp;
            stat(fLoad.str().c_str(), &statTmp);

            // Check if it is a file or directory
            if( !S_ISDIR(statTmp.st_mode) )
            {
                cout << "Loading file " << fLoad.str() << endl;

                LoadSingleGraph(string(fLoad.str()));
                fLoad.str("");
                iCounter++;
            }
            else
            {
                cout << "Loading directory " << fLoad.str() << endl;
                iCounter += LoadAllGraphs(string(fLoad.str()));
                fLoad.str("");
            }
        }
        if( iCounter != -1 && iCounter >= iLimit )
            return iCounter;
    }

    closedir(dp);

    return iCounter;
}
//----------------------------------------------------------------
void GraphHandler::ListMetaVertices() const
{
    for(pair<MetaIter_t, MetaIter_t> p = vertices(m_MetaGraph); p.first!=p.second; ++p.first)
        cout << "Vertex: " << m_MetaGraph[*p.first].sName << ", Freq: " << m_MetaGraph[*p.first].iFreq << endl;
}
//----------------------------------------------------------------
void GraphHandler::SaveMetaGraph(string sFilename)
{

   // if( m_bMetaFiltered )
   //     SaveMetaGraphInternal(sFilename, *m_pFiltMGraph);
   // else
        SaveMetaGraphInternal(sFilename, m_MetaGraph);

}
//----------------------------------------------------------------
template<typename GeneralGraph>
void GraphHandler::SaveMetaGraphInternal(string sFilename, const GeneralGraph& G) const
{

    // Create the vertex and edge writers
    MetaVertex_writer<GeneralGraph> vertexwriter(G);
    MetaEdge_writer<GeneralGraph> edgewriter(G);

    // Save the graph to the file
    std::ofstream ofs( sFilename.c_str() );
    boost::write_graphviz(ofs, G, vertexwriter, edgewriter );
}
//----------------------------------------------------------------
void GraphHandler::SaveGraph(string sFilename, const Graph& G, vector<const Graph*> pOverlay) const
{


    // Save the graph to the file
    std::ofstream ofs( sFilename.c_str() );
    if( pOverlay.size() == 0)
    {
        // Create the vertex and edge writers
        Vertex_writer vertexwriter(G);

        boost::write_graphviz(ofs, G, vertexwriter );
    }
    else
    {
        // Create the vertex and edge writers
        Vertex_writer vertexwriter(G, pOverlay);
        Edge_writer<Graph> edgewriter(G, pOverlay);

        boost::write_graphviz(ofs, G, vertexwriter, edgewriter );

    }
}
//----------------------------------------------------------------
void GraphHandler::SaveGraph(string sFilename, const Graph& G, const Graph* pOverlay) const
{
    vector<const Graph*> kTmp;
    if( pOverlay )
        kTmp.push_back(pOverlay);
    SaveGraph(sFilename, G, kTmp);
}
//----------------------------------------------------------------
void GraphHandler::FilterMetaGraph(int iVertexCutoff, int iEdgeCutoff)
{
/*
    if( !m_bMetaFiltered )
    {
        // Create two filter objects
        m_Vertexfilter = vertex_freq_filter<MetaGraph>(m_MetaGraph, iVertexCutoff);
        m_Edgefilter = edge_freq_filter<MetaGraph>(m_MetaGraph, iEdgeCutoff);

        if( m_pFiltMGraph )
            delete m_pFiltMGraph;

        // Create the filtered graph with the objects defined above
        m_pFiltMGraph = new FiltMGraph(m_MetaGraph, m_Edgefilter, m_Vertexfilter);

        m_bMetaFiltered = true;
    }
    */
}
void GraphHandler::Blacklist(const set<string>& sBlacklist)
{
    BOOST_FOREACH( Graph& G, m_Graphs )
        G = BlacklistInternal(G, sBlacklist);
}
//----------------------------------------------------------------
GraphHandler::Graph GraphHandler::BlacklistInternal
            (const GraphHandler::Graph& G, const set<string>& sBlacklist)
{
    // We shall constuct a temporary new graph and then
    // overwrite the old
    Graph gTmp;



    Iter_t vi, vi_end, next;

    // Maintain a mapping from the old vertex descriptor -> new descriptor
    map< Vertex, Vertex*> tmpMapping;

    bool bBlacklisted = false;
    for (tie(vi, vi_end) = vertices(G); vi != vi_end; ++vi)
    {
        bBlacklisted = false;
        tmpMapping[*vi] = 0;

        //bBlacklisted =
        //    (find(sBlacklist.begin(), sBlacklist.end(), G[*vi].sName) != sBlacklist.end());
        BOOST_FOREACH( const string& s, sBlacklist )
        {

            if( m_IDTypes[get(vertex_name, G, *vi)].compare(s) == 0 )
                bBlacklisted = true;
        }

        if( !bBlacklisted )
        {
            Vertex* ver = new Vertex(add_vertex(gTmp));
            //gTmp[*ver] = G[*vi];
            int iID = get(vertex_name, G, *vi);//m_TypesID.find(G[*vi].sName)->second;

            put(vertex_name, gTmp, *ver, iID);

            tmpMapping[*vi] = ver;
        }
    }


    graph_traits<Graph>::edge_iterator firstE, lastE;
    // Go through all the edges of the old graph and only add those which have verticies
    // which are existing in the transformed graph.

    // Also, for all those edges (u,v) where u or v is one of the blacklisted vertices
    // make it so that its edges
    for (tie(firstE, lastE) = edges(G); firstE != lastE; ++firstE)
    {
            if( tmpMapping[source(*firstE, G)] != 0 && tmpMapping[target(*firstE, G)] != 0  )
            {

                pair<Edge, bool> kRet2 = edge(*tmpMapping[source(*firstE, G)], *tmpMapping[target(*firstE, G)], gTmp);

                if( !kRet2.second )
                    std::pair<Edge, bool> kRet =
                    add_edge(*tmpMapping[source(*firstE, G)], *tmpMapping[target(*firstE, G)], gTmp);
            }
            else
            {
                // Either u or v is blacklisted, (or both).
                // If both are blacklisted, we want to connect all the valid adjacent vertices
                // to each other

                // Assume only one is blacklisted
                Vertex vBlackList;
                Vertex vOther;
                if( tmpMapping[source(*firstE, G)] == 0 )
                {
                    vBlackList = source(*firstE, G);
                    vOther = target(*firstE, G);
                }
                else
                {
                    vBlackList = target(*firstE, G);
                    vOther = source(*firstE, G);
                }

                if( tmpMapping[vOther] != 0 )
                {
                // Go through all the blacklisted ones adjacent vertices
                // We want to connect these to the vOther vertex (if they are not already connected)
                graph_traits<Graph>::adjacency_iterator v1, v2;
                bool bFound = false;
                for( tie(v1, v2) = adjacent_vertices(vBlackList, G); (v1 != v2); ++v1 )
                {
                    // Check if this vertex *v1 is already connected to VOther
                    // If not, add the connection
                    pair<Edge, bool> kRet = edge( *v1, vOther, G);

                    // Add if it doesn't exist
                    if( !kRet.second && *v1 != vOther )
                        kRet = add_edge(*tmpMapping[*v1], *tmpMapping[vOther], gTmp);
                }
                }
            }
    }

    //Graph kFinalGraph;

    for(map<Vertex, Vertex*>::const_iterator i = tmpMapping.begin(); i != tmpMapping.end(); ++i)
        delete i->second;

    tmpMapping.clear();

    // 3. We must also remove those vertices which have no more edges
    //---------------------------------------------
    RemoveIsolatedVertices(&gTmp);
    //---------------------------------------------
    // End of step 3

    // Copy the graph name as well
    set_property(gTmp, boost::graph_name, get_property(G,boost::graph_name));

    return gTmp;
}
//----------------------------------------------------------------
// Removes those vertices which have no edges
void GraphHandler::RemoveIsolatedVertices()
{
    BOOST_FOREACH(Graph&G, m_Graphs)
            RemoveIsolatedVertices(&G);
}

// Removes those vertices which have no edges
void GraphHandler::RemoveIsolatedVertices(Graph* pG)
{
    Graph kFinalGraph;
    Iter_t vi, vi_end;
    graph_traits<Graph>::edge_iterator firstE, lastE;

    Graph& G = *pG;

    // Maintain a mapping from the old vertex descriptor -> new descriptor
    map< Vertex, Vertex*> tmpMapping;
//out_edges(v, g)
    for (tie(vi, vi_end) = vertices(G); vi != vi_end; ++vi)
    {
        tmpMapping[*vi] = 0;
        bool bC; bC = true;

        if( out_degree(*vi, G) == 1 )
        {
            boost::graph_traits<Graph>::out_edge_iterator e, e_end;

            tie(e, e_end) = out_edges(*vi, G);
            //get(property, g, x)
            if( get(vertex_name, G, target(*e, G)) == get(vertex_name, G, *vi))
                bC = false;
        }


        // Only add if this vertex has some edges, and if that connection is not to itself only.
        if( out_degree(*vi, G) != 0 && bC )
                //&& !(   out_degree(*vi, G) == 1 &&  (target(*out_edges(*vi, G).first, G)==(*vi))   ))
        {
            // Add it
            Vertex* kV = new Vertex(add_vertex(kFinalGraph));
            //kFinalGraph[*kV] = G[*vi];
            int iID = get(vertex_name, G, *vi);

            put(vertex_name, kFinalGraph, *kV, iID);

            tmpMapping[*vi] = kV;
        }
    }


    // Transfer all edges
    for (tie(firstE, lastE) = edges(G); firstE != lastE; ++firstE)
        if( tmpMapping[source(*firstE, G)] != 0 && tmpMapping[target(*firstE, G)] != 0  )
            add_edge(*tmpMapping[source(*firstE, G)], *tmpMapping[target(*firstE, G)], kFinalGraph);

    for(map<Vertex, Vertex*>::const_iterator i = tmpMapping.begin(); i != tmpMapping.end(); ++i)
        delete i->second;

    set_property(kFinalGraph, boost::graph_name, get_property(G,boost::graph_name));
    G = kFinalGraph;
}
//----------------------------------------------------------------
void GraphHandler::FilterTransformMetaGraph(int iVertexCutoff, int iEdgeCutoff)
{
    MetaIter_t vi, vi_end, next;

    // We shall constuct a temporary new graph and then
    // overwrite the old m_MetaGraph
    MetaGraph gTmp;

    // Maintain a mapping from the old vertex descriptor -> new descriptor
    map< MetaVertex, MetaVertex*> tmpMapping;

    // We will blacklist all those vertices which are below the threshold
    // from the set of all graphs.
    set<string> kToBlacklist;

    for (tie(vi, vi_end) = vertices(m_MetaGraph); vi != vi_end; ++vi)
    {
        tmpMapping[*vi] = 0;



        if( m_MetaGraph[*vi].iFreq > iVertexCutoff )
        {
            MetaVertex* ver = new MetaVertex(add_vertex(gTmp));
            gTmp[*ver] = m_MetaGraph[*vi];

            tmpMapping[*vi] = ver;
        }
        else
            kToBlacklist.insert(m_MetaGraph[*vi].sName);
    }

    graph_traits<MetaGraph>::edge_iterator firstE, lastE;
    // Gh through all the edges of the old graph and only add those which have verticies
    // which are existing in the transformed graph.
    for (tie(firstE, lastE) = edges(m_MetaGraph); firstE != lastE; ++firstE)
        if( m_MetaGraph[*firstE].iFreq > iEdgeCutoff )
            if( tmpMapping[source(*firstE, m_MetaGraph)] != 0 && tmpMapping[target(*firstE, m_MetaGraph)] != 0  )
            {
                std::pair<MetaEdge, bool> kRet =
                add_edge(*tmpMapping[source(*firstE, m_MetaGraph)], *tmpMapping[target(*firstE, m_MetaGraph)], gTmp);

                gTmp[kRet.first] = m_MetaGraph[*firstE];
            }

    for(map<MetaVertex, MetaVertex*>::const_iterator i = tmpMapping.begin(); i != tmpMapping.end(); ++i)
        delete i->second;

    m_MetaGraph = gTmp;

    Blacklist(kToBlacklist);
   // MergeVertices(kToBlacklist, "Others");

    // Go through all the graphs and remove those which have 0 vertices.
    RemoveEmptyGraphs();
}

void GraphHandler::RemoveEmptyGraphs()
{
    vector<Graph> kTmpV;
    int iX = 0;
    BOOST_FOREACH( const Graph& G, m_Graphs )
    {
        if( num_vertices(G) != 0 )
        {
            kTmpV.push_back(G);
        }
        iX++;
    }
    m_Graphs = kTmpV;
}
void GraphHandler::MergeVertices(const set<string>& ToMerge, string sNewName)
{
    // 1. Create a new vertex sNewName in each graph that has either sName1 or sName2 in it.
    // 2. Traverse all the edges of sName1 and sName2 and connect those nodes to sNewName
    // 3. Remove sName1 and sName2

    BOOST_FOREACH( Graph& G, m_Graphs )
        MergeVerticesInternal(ToMerge, G, sNewName);

    // Do the same for the metagraph, except maintain the internal properties as well.
    //if( m_bMetaFiltered )
    //    MergeVerticesInternal(ToMerge, *m_pFiltMGraph, sNewName);
    //else
        MergeVerticesInternalMeta(ToMerge, m_MetaGraph, sNewName);

}

void GraphHandler::MergeVerticesInternal(const set<string>& ToMerge,
                Graph& G, string sNewName)
{
    bool bAdded = false;

    typedef graph_traits<Graph>::vertex_iterator GnrItr;

    graph_traits<Graph>::vertex_descriptor pNew;

    pair<bool, GnrItr> kRet;

        BOOST_FOREACH( const string& sName, ToMerge )
        {
            if( (kRet = FindVertex(sName, G)).first )
            {
                // This vector(sName) exists in this specific graph
                if( !bAdded )
                {
                    pNew = add_vertex(G);

                    //G[pNew].sName = sNewName;
                    //get(vertex_name, G, *vi)

                    m_TypesID[sNewName] = m_iIDCounter;
                    m_IDTypes[m_iIDCounter] = sNewName;
                    m_iIDCounter++;

                    put(vertex_name, G, pNew, m_TypesID[sNewName]);


                    bAdded = true;

                    // Invalidated?
                    kRet = FindVertex(sName, G);
                }

                typedef graph_traits<Graph>::out_edge_iterator Edg_It;
                Edg_It first, last;


                // Add all the edges to the new vertex
                for (tie(first, last) = out_edges(*kRet.second, G); first != last; ++first)
                {

                    pair<Edge, bool> kRet = edge( pNew, target(*first, G), G);

                    // Add if it doesn't exist
                    if( !kRet.second )
                        kRet = add_edge(pNew, target(*first, G), G);

                }
                clear_vertex(*kRet.second,  G); // Remove all edges from the vertex
                remove_vertex(*kRet.second, G); // Remove the vertex from the graph

                pNew = *FindVertex(sNewName, G).second;
            }
        }
}
// TODO: This specialization is needed at the moment since MetaGraph
// has bundled properties which the normal graphs do not have (and these needs to be preserved)
// when merging. However, alot of the code is very similar, so perhaps there are other methods of
// code reuse that could be used in this case. This will have to do for now.
template<typename GeneralGraph>
void GraphHandler::MergeVerticesInternalMeta(const set<string>& ToMerge,
                GeneralGraph& G, string sNewName)
{
    /*
    bool bAdded = false;

    typedef typename graph_traits<GeneralGraph>::vertex_iterator GnrItr;

    typename graph_traits<GeneralGraph>::vertex_descriptor pNew;

    pair<bool, GnrItr> kRet;

        BOOST_FOREACH( const string& sName, ToMerge )
        {
            if( (kRet = FindVertex(sName, G)).first )
            {
                // This vector(sName) exists in this specific graph
                if( !bAdded )
                {
                    pNew = add_vertex(G);

                    //put(vertex_name, G, pNew, sNewName);
                    G[pNew].sName = sNewName;
                    G[pNew].iFreq += G[*kRet.second].iFreq;
                    bAdded = true;
                }
                else
                {
                    G[pNew].iFreq += G[*kRet.second].iFreq;
                }

                typedef typename graph_traits<GeneralGraph>::out_edge_iterator Edg_It;
                Edg_It first, last;

                // Add all the edges to the new vertex
                for (tie(first, last) = out_edges(*kRet.second, G); first != last; ++first)
                {


                    pair<MetaEdge, bool> kRet = edge( pNew, target(*first, G), m_MetaGraph);

                    // Add if it doesn't exist
                    if( !kRet.second )
                        kRet = add_edge(pNew, target(*first, G), G);

                    G[kRet.first].iFreq += G[*first].iFreq;
                }
                clear_vertex(*kRet.second,  G); // Remove all edges from the vertex
                remove_vertex(*kRet.second, G); // Remove the vertex from the graph

                // Since all vertex_descriptors are invalid
                pNew = *FindVertex(sNewName, G).second;

                //cout << "Removed one vertex" << endl;
                //ListMetaVertices();
            }
        }
        */
}
//----------------------------------------------------------------
GraphHandler::Graph& GraphHandler::GetGraph(int iID)
{
    return m_Graphs[iID];
}
//----------------------------------------------------------------
int GraphHandler::GetNumGraphs() const
{
    return m_Graphs.size();
}
//----------------------------------------------------------------
// Makes sure that no graphs are isomorphic to each other
void GraphHandler::RemoveDuplicates(std::vector<Graph>* pVec)
{
    if( pVec == 0 )
        pVec = &m_Graphs;

    vector<Graph> kNewGraphs;
    vector<bool> bChecked; // Keep track on which entries we've already checked
    for( unsigned int i = 0; i < pVec->size(); i++ )
        bChecked.push_back(false);



    for( unsigned int i = 0; i < pVec->size(); i++ )
    {
        // If we haven't checked this coordinate yet
        if( !bChecked[i])
        {
            bChecked[i] = true;
            kNewGraphs.push_back((*pVec)[i]);

            for( unsigned int j = i + 1; j < pVec->size(); j++ )
            {
                // They are equivalent
                //if( !bChecked[j] && IsoGraphEditDistance::GraphEditDistance((*pVec)[i], (*pVec)[j]) == 0 )
                if( !bChecked[j] && EditDistance((*pVec)[i], (*pVec)[j]) == 0 )
                {
                    //cout << "Graph " << i << " and graph " << j << " are eqv." << endl;
                    bChecked[j] = true;
                }
            }
        }
    }
    //cout << "Old size " << m_Graphs.size() << endl;
    //cout << "New size " << kNewGraphs.size() << endl;
    pVec->clear();
    *pVec = kNewGraphs;
}
//----------------------------------------------------------------
// Implementation of Algorithm 4 from "Efficient Algorithms on Trees and Graphs with Unique
//Node Labels"
// See also "On Graphs with Unique Node Labels" Bunke
int GraphHandler::EditDistance(const Graph& G1, const Graph& G2) const
{
    int iDistance = 0;

    // 1. Sort V1 and V2 by node label
    Iter_t G1begin, G1end; tie(G1begin, G1end) = vertices(G1);
    Iter_t G2begin, G2end; tie(G2begin, G2end) = vertices(G2);

    vector<Vertex> V1, V2;

    V1.resize(num_vertices(G1));
    V2.resize(num_vertices(G2));

    copy(G1begin, G1end, V1.begin());
    copy(G2begin, G2end, V2.begin());


    Sorter kSortG1(G1);
    Sorter kSortG2(G2);

    sort(V1.begin(), V1.end(), kSortG1);
    sort(V2.begin(), V2.end(), kSortG2);



    // 2.
    while( V1.size() != 0 && V2.size() != 0 )
    {
        // v1, v2 first vertices of V1, V2 respectivly
        // if alpha1(v1) < alpha2(v2)
        int a1 = get(vertex_name, G1, *V1.begin());
        int a2 = get(vertex_name, G2, *V2.begin());
        if( a1 < a2 )
        {
            iDistance++;
           // cout << "Delete node" << get(vertex_name, G1, *V1.begin()) << endl;
            V1.erase(V1.begin());
        }
        else if( a1 > a2 )
        {
           // cout << "Insert node" << get(vertex_name, G2, *V2.begin()) << endl;
            iDistance++;
            V2.erase(V2.begin());
        }
        else
        {
            V1.erase(V1.begin());
            V2.erase(V2.begin());
        }
    }
    if( V1.size() != V2.size() )
    {
        iDistance += V1.size() > V2.size() ? V1.size() - V2.size() : V2.size() - V1.size();
    }


    // 3. Sort E1 and E2 by target node label
    // Sort E1 and E2 by source node label
    graph_traits<Graph>::edge_iterator G1Ebegin, G1Eend, G2Ebegin, G2Eend;
    tie(G1Ebegin, G1Eend) = edges(G1); tie(G2Ebegin, G2Eend) = edges(G2);



    vector<Edge> E1, E2;
    E1.resize(num_edges(G1));
    E2.resize(num_edges(G2));

    copy(G1Ebegin, G1Eend, E1.begin());
    copy(G2Ebegin, G2Eend, E2.begin());

    EdgeSorter eSortG1(G1);
    EdgeSorter eSortG2(G2);
    sort(E1.begin(), E1.end(), eSortG1);
    sort(E2.begin(), E2.end(), eSortG2);

    eSortG1.bTargetNode = false; eSortG2.bTargetNode = false;
    sort(E1.begin(), E1.end(), eSortG1);
    sort(E2.begin(), E2.end(), eSortG2);

    //4.
    while( E1.size() != 0 && E2.size() != 0 )
    {
        // Let (v1, w1) and (v2, w2) be the first elements of E1, E2
        // source(*E1.begin()) == v1
        // target(*E1.begin()) == w1

        // source(*E2.begin()) == v2
        // target(*E2.begin()) == w2

        int a1 = get(vertex_name, G1, source(*E1.begin(), G1)); // alpha1(v1)
        int a2 = get(vertex_name, G1, source(*E2.begin(), G2)); // alpha2(v2)

        int w1 = get(vertex_name, G1, target(*E1.begin(), G1)); // alpha1(w1)
        int w2 = get(vertex_name, G1, target(*E2.begin(), G2)); // alpha2(w2)


        // If alpha1(v1) < alpha2(v2)
        if( a1 < a2 || (a1==a2 && w1 < w2) )
        {
            //cout << "Delete edge " << a1 << " to " << w1;
            iDistance++;
            E1.erase(E1.begin());
        }
        else if( a1 > a2 || ( a1==a2 && w1 > w2 ) )
        {
           // cout << "Insert edge " << a2 << " to " << w2;
            iDistance++;
            E2.erase(E2.begin());
        }
        else
        {
            E1.erase(E1.begin());
            E2.erase(E2.begin());
        }
    }
    if( E1.size() != E2.size() )
    {
        iDistance += E1.size() > E2.size() ? E1.size() - E2.size() : E2.size() - E1.size();
    }

    return iDistance;
}
//----------------------------------------------------------------
pair<bool, GraphHandler::Vertex> GraphHandler::FindVertex(unsigned int iName, const Graph& G) const
{
    GraphHandler::Iter_t vi, vi_end;
    // For each vertex
    for (tie(vi, vi_end) = vertices(G); vi != vi_end; ++vi)
    {
        if( get(vertex_name, G, *vi) == iName )
            return make_pair(true, *vi);
    }

    return make_pair(false, Vertex());
}
//----------------------------------------------------------------
// Check if G1 is contained in G2
bool GraphHandler::IsSubgraphIsomorph(const Graph& G1, const Graph& G2) const
{
    // Pointer to the smallest graph
    const Graph* pSG;
    const Graph* pLG;

/*
    if( num_vertices(G1) <  num_vertices(G2) )
    {
        pSG = &G1;
        pLG = &G2;
    }
    else
    {
        pSG = &G2;
        pLG = &G1;
    }*/

    pSG = &G1;
    pLG = &G2;

    map<unsigned int, Vertex> kMap;

    // Check that all the vertices of the smallest graph and all the edges
    // are present in the bigger graph


    GraphHandler::Iter_t vi, vi_end;
    GraphHandler::Iter_t vi1, vi_end1;

    for (tie(vi, vi_end) = vertices(*pSG); vi != vi_end; ++vi)
    {
        bool bFound = false;
        // Check if we can find this vertex in the larger graph
        for (tie(vi1, vi_end1) = vertices(*pLG); vi1 != vi_end1; ++vi1)
        {
            if( get(vertex_name, *pSG, *vi) == get(vertex_name, *pLG, *vi1) )
            {
                kMap[get(vertex_name, *pLG, *vi1)] = *vi1;
                bFound = true;
            }
        }
        if(!bFound)
            return false;
    }

    for (tie(vi, vi_end) = vertices(*pSG); vi != vi_end; ++vi)
    {
        for (tie(vi1, vi_end1) = vertices(*pSG); vi1 != vi_end1; ++vi1)
        {
            pair<Edge, bool> kRet1 = edge( *vi, *vi1, *pSG);
            pair<Edge, bool> kRet2 = edge( kMap[get(vertex_name, *pSG, *vi)], kMap[get(vertex_name, *pSG, *vi1)], *pLG);

            // If the edge exist in the smaller graph, it must also exist in the larger graph
            if( kRet1.second )
                if( !kRet2.second )
                    return false;
        }
    }

    return true;
}

// Based on Algorithm 1 in Efficient Algorithms on Trees and Graphs with Unique
// Node Labels
GraphHandler::Graph GraphHandler::LargestCommonSubGraph(Graph& G1, Graph& G2)
{
    Graph gMCS;

    GraphHandler::Iter_t vi, vi_end;
    GraphHandler::Iter_t vi1, vi_end1;

    map<Vertex, Vertex> G1Map; // Maps the new vertex into the corresponding "old" from G1
    map<Vertex, Vertex> G2Map; // Maps the new vertex into the corresponding "old" from G1


    // For each vertex in G1
    for (tie(vi, vi_end) = vertices(G1); vi != vi_end; ++vi)
    {
        for (tie(vi1, vi_end1) = vertices(G2); vi1 != vi_end1; ++vi1)
        {
            if( get(vertex_name, G1, *vi) == get(vertex_name, G2, *vi1) )
            {
                Vertex vNew;
                // Add it
                vNew = add_vertex(gMCS);
                //gMCS[vNew].sName = G1[*vi].sName;
                put(vertex_name, gMCS, vNew, get(vertex_name, G1, *vi));

                G1Map[vNew] = *vi;
                G2Map[vNew] = *vi1;
            }
        }
    }

    for (tie(vi, vi_end) = vertices(gMCS); vi != vi_end; ++vi)
    {
        // let v1 be the node of G1 with α1 (v1 ) = α(v)
        // In our case, v1=G1Map[*vi], v2=G2Map[*vi]
        for (tie(vi1, vi_end1) = vertices(gMCS); vi1 != vi_end1; ++vi1)
        {
           // let w1 be the node of G1 with α1 (w1 ) = α(w)
            // w1=G1Map[*vi1], w2=G2Map[*vi1]
            // if (v1 , w1 ) ∈ E1 , (v2 , w2 ) ∈ E2 and β1 (v1 , w1 ) = β2 (v2 , w2 ) then

            //Vertex vSource = kMap[get(vertex_name, G, source(*firstE, G))];
            //Vertex vTarget = kMap[get(vertex_name, G, target(*firstE, G))];

            pair<Edge, bool> kRet1 = edge( G1Map[*vi], G1Map[*vi1], G1);
            pair<Edge, bool> kRet2 = edge( G2Map[*vi], G2Map[*vi1], G2);

            if( kRet1.second && kRet2.second )
            {
                //E := E ∪ {(v, w)}
                add_edge(*vi, *vi1, gMCS);
            }
        }
    }
    RemoveIsolatedVertices(&gMCS);
    return gMCS;

    /*
    int iDistance = 0;

    Graph gMCS;

    // 1. Sort V1 and V2 by node label
    Iter_t G1begin, G1end; tie(G1begin, G1end) = vertices(G1);
    Iter_t G2begin, G2end; tie(G2begin, G2end) = vertices(G2);

    vector<Vertex> V1, V2;

    V1.resize(num_vertices(G1));
    V2.resize(num_vertices(G2));

    copy(G1begin, G1end, V1.begin());
    copy(G2begin, G2end, V2.begin());

    Sorter kSortG1(G1);
    Sorter kSortG2(G2);

    sort(V1.begin(), V1.end(), kSortG1);
    sort(V2.begin(), V2.end(), kSortG2);

    Vertex vNew;

    map<Vertex, Vertex> gamma;

    // 2.
    while( V1.size() != 0 && V2.size() != 0 )
    {
        // v1, v2 first vertices of V1, V2 respectivly
        // if alpha1(v1) < alpha2(v2)
        int a1 = get(vertex_name, G1, *V1.begin());
        int a2 = get(vertex_name, G2, *V2.begin());
        if( a1 == a2 )
        {
            iDistance++;
           // cout << "Delete node" << get(vertex_name, G1, *V1.begin()) << endl;
            V1.erase(V1.begin());
        }
        else if( a1 > a2 )
        {
           // cout << "Insert node" << get(vertex_name, G2, *V2.begin()) << endl;
            iDistance++;
            V2.erase(V2.begin());
        }
        else
        {
            vNew = add_vertex(gMCS);
            gMCS[vNew].sName = G1[*V1.begin()].sName;
            put(vertex_name, gMCS, vNew, get(vertex_name, G1, *V1.begin()));

            gamma[*V1.begin()] = vNew; // Record it

            V1.erase(V1.begin());
            V2.erase(V2.begin());
        }
    }
    if( V1.size() != V2.size() )
    {
        iDistance += V1.size() > V2.size() ? V1.size() - V2.size() : V2.size() - V1.size();
    }


    // 3. Sort E1 and E2 by target node label
    // Sort E1 and E2 by source node label
    graph_traits<Graph>::edge_iterator G1Ebegin, G1Eend, G2Ebegin, G2Eend;
    tie(G1Ebegin, G1Eend) = edges(G1); tie(G2Ebegin, G2Eend) = edges(G2);



    vector<Edge> E1, E2;
    E1.resize(num_edges(G1));
    E2.resize(num_edges(G2));

    copy(G1Ebegin, G1Eend, E1.begin());
    copy(G2Ebegin, G2Eend, E2.begin());

    EdgeSorter eSortG1(G1);
    EdgeSorter eSortG2(G2);
    sort(E1.begin(), E1.end(), eSortG1);
    sort(E2.begin(), E2.end(), eSortG2);

    eSortG1.bTargetNode = false; eSortG2.bTargetNode = false;
    sort(E1.begin(), E1.end(), eSortG1);
    sort(E2.begin(), E2.end(), eSortG2);

    //4.
    while( E1.size() != 0 && E2.size() != 0 )
    {
        // Let (v1, w1) and (v2, w2) be the first elements of E1, E2
        // source(*E1.begin()) == v1
        // target(*E1.begin()) == w1

        // source(*E2.begin()) == v2
        // target(*E2.begin()) == w2

        int a1 = get(vertex_name, G1, source(*E1.begin(), G1)); // alpha1(v1)
        int a2 = get(vertex_name, G1, source(*E2.begin(), G2)); // alpha2(v2)

        int w1 = get(vertex_name, G1, target(*E1.begin(), G1)); // alpha1(w1)
        int w2 = get(vertex_name, G1, target(*E2.begin(), G2)); // alpha2(w2)

        //if()

        // If alpha1(v1) < alpha2(v2)
        if( a1 < a2 || (a1==a2 && w1 < w2) )
        {
            //cout << "Delete edge " << a1 << " to " << w1;
            iDistance++;
            E1.erase(E1.begin());
        }
        else if( a1 > a2 || ( a1==a2 && w1 > w2 ) )
        {
           // cout << "Insert edge " << a2 << " to " << w2;
            iDistance++;
            E2.erase(E2.begin());
        }
        else
        {
            //gamma
            add_edge(gamma[source(*E1.begin(), G1)], gamma[target(*E1.begin(), G1)], gMCS);
            E1.erase(E1.begin());
            E2.erase(E2.begin());
        }
    }
    if( E1.size() != E2.size() )
    {
        iDistance += E1.size() > E2.size() ? E1.size() - E2.size() : E2.size() - E1.size();
    }

    //RemoveIsolatedVertices(&gMCS);
    return gMCS;*/
}
//----------------------------------------------------------------
void GraphHandler::SaveGSpanData(string sFilename, const vector<Graph>& kVec)

{

    ofstream out(sFilename.c_str(), ios_base::out);

    for(int i=0; i<kVec.size(); i++)
    {
        out << "t # " << i << endl;
        const Graph& G = kVec[i];


        GraphHandler::Iter_t vi, vi_end;
        int j = 0;
        // For each vertex
        for (tie(vi, vi_end) = vertices(G); vi != vi_end; ++vi)
        {
            out << "v " << j << " " << get(vertex_name, G, *vi) << endl;
            j++;
        }

        graph_traits<GraphHandler::Graph>::edge_iterator firstE, lastE;
        // For each edge
        for (tie(firstE, lastE) = edges(G); firstE != lastE; ++firstE)
        {
            out << "e " << boost::source(*firstE, G) << " "
                         << boost::target(*firstE, G) << " 1" << endl;

        }
    }



    out.close();
}

vector<GSpanStruct> GraphHandler::LoadGSpanData(string sFilename)
{

    ifstream in(sFilename.c_str());

    vector<GSpanStruct> kVecTmp;


    string sTmp;
    while( !in.eof() )
    {

        GSpanStruct gTmp;
        Graph& G = gTmp.G;

        std::vector<std::string> strs;

        getline(in, sTmp, '\n');
        boost::split(strs, sTmp, boost::is_any_of("\t "));


        while( strs.size() != 0 && strs[0].size() != 0 )
        {

            if( strs[0].compare("t") == 0 )
            {
                //t # 0 * 272
                gTmp.iSupport = lexical_cast<int>(strs[4]);

            }

            // Vertex found
            if( strs[0].compare("v") == 0 )
            {
                Vertex vNew = add_vertex(G);
                int iID = lexical_cast<int>(strs[2]);

                //G[vNew].sName = m_IDTypes[iID];
                put(vertex_name, G, vNew, iID);

            }

            if( strs[0].compare("e") == 0 )
            {
                // Add an edge
                add_edge(lexical_cast<Vertex>(strs[1]), lexical_cast<Vertex>(strs[2]), G);
            }

            if( strs[0].compare("x") == 0 )
            {
                for(int l = 1; l<strs.size()-1; l++)
                {
                    gTmp.kGraphIDs.push_back(lexical_cast<int>(strs[l]));
                }
            }

            getline(in, sTmp, '\n');
            boost::split(strs, sTmp, boost::is_any_of("\t "));

        }
        gTmp.G = G;
        kVecTmp.push_back(gTmp);
    }

    return kVecTmp;
}


vector<GSpanStruct> GraphHandler::GSpan(const vector<Graph>& kVec, double dFreq)
{
    SaveGSpanData("tmp", kVec);

    stringstream sTmp;

    sTmp << "../thirdparty/gSpan6/gSpan -f tmp -s " << dFreq << " -o -i > output.txt";

    //system ("../thirdparty/gSpan6/gSpan -f tmp -s 0.1 -o -i");
    system (sTmp.str().c_str());
    return LoadGSpanData("tmp.fp");
}

vector<GSpanStruct> GraphHandler::GSpan(const vector<Graph*>& kVec, double dFreq)
{
    vector<Graph> kVecTmp;

    for(int i=0; i<kVec.size(); i++)
    {
        kVecTmp.push_back(*(kVec[i]));
    }
    return GSpan(kVecTmp, dFreq);
}

// Returns true if we can find some resulting subgraph that has minSup as support (the number of times
// that subgraph is occuring among the set of Graphs), and has a size >= minSize
pair<bool, GraphHandler::Graph> GraphHandler::GSpanStructure(const vector<Graph*>& kVec, int minSup, int minSize)
{
    vector<Graph> kVecTmp;

    for(int i=0; i<kVec.size(); i++)
    {
        kVecTmp.push_back(*(kVec[i]));
    }

    const vector<GSpanStruct>& kRes = GSpan(kVecTmp, 1);

    for(int i=0; i<kRes.size(); i++)
    {
        if( kRes[i].iSupport >= minSup && num_vertices(kRes[i].G) >= minSize )
        {
            //return pair<bool, Graph>(true, kRes[i].G);
            int jb = i;
            int iMax = num_vertices(kRes[i].G);
            for(int j=0; j<kRes.size(); j++)
            {
                if( num_vertices(kRes[j].G) > iMax )
                    jb = j;
            }
            return pair<bool, Graph>(true, kRes[jb].G);
        }
    }
    return pair<bool, Graph>(false, Graph());
}

// Implementation of Algorithm 1 from Online Structural Graph Clustering Using Frequent subgraph mining
vector< vector<GraphHandler::Graph*> > GraphHandler::StructureClustering(double dTheta, int iMinGraphSize)
{
    vector< vector<Graph*> > clusters;

    // Used to keep track of the graph in the cluster with the maximum/minimum num of vertices
    vector< int > ClusterMax;
    vector< int > ClusterMin;

    kCommon.clear();

    int iMinSup;

    //for(unsigned int meta=0; meta<2; meta++)
   // {
 //   iMinGraphSize++;
    // Loop over all graphs
    for(unsigned int j=0; j<m_Graphs.size(); j++)
    {
        cout << "Graph "<< j << endl;
        bool bHasCluster = false;
        int iSizeQuery = num_vertices(m_Graphs[j]); // |x_m+1|

        if( iSizeQuery >= iMinGraphSize )
        {
            // Compare graph against existing clusters
            for(unsigned int c=0; c < clusters.size(); c++)
            {
                // max(size(graph[j]), size(c.max))
                int dMinSize = (int)(dTheta * (double)max(iSizeQuery, ClusterMax[c]));

                if( ClusterMax.size()-1 < c )
                {
                    ClusterMax.push_back(0);
                    ClusterMin.push_back(0);
                }

                // x_m+1 is the new graph (query graph)
                // x_max is the largest graph in the cluster
                // x_min is the smallest graph in the cluster
                // (3):  |x_m+1| > |x_max| ^ minSize > |x_min|
                // (4):  |x_m+1| < |x_min| ^ minSize > |x_m+1|
                bool b3 = iSizeQuery > ClusterMax[c] && dMinSize > ClusterMin[c];
                bool b4 = iSizeQuery < ClusterMin[c] && dMinSize > iSizeQuery;

                // Check (3) || (4)
                if( b3 || b4 )
                {
                    continue;
                }
                else
                {
                    // We require the support to be the whole set (cluster union query graph)
                    // That is, we require it exists in all the graphs.
                    iMinSup = clusters[c].size() + 1;

                    vector<Graph*> kVec;
                    kVec.resize(clusters[c].size());
                    copy(clusters[c].begin(), clusters[c].end(), kVec.begin());
                    kVec.push_back(&m_Graphs[j]);

                    // add graph[j] to cluster c if gSpan finds at least one
                    // // common subgraph that meets the minSize threshold
                    // if gSpan (graph[j] ∪ c.graphs, minSup, minSize) then
                    pair<bool, Graph> kPTmp = GSpanStructure(kVec, iMinSup, iMinGraphSize);
                    if( kPTmp.first )
                    {

                        kCommon[c] = kPTmp.second;
                        // Add to cluster
                        if( find(clusters[c].begin(), clusters[c].end(), &m_Graphs[j]) == clusters[c].end() )
                        {
                            clusters[c].push_back(&m_Graphs[j]);

                        // Update max/min as well
                        if( iSizeQuery > ClusterMax[c] )
                            ClusterMax[c] = iSizeQuery;
                        if( iSizeQuery < ClusterMin[c] )
                            ClusterMin[c] = iSizeQuery;
                        bHasCluster = true;
                        }
                    }
                }

            } // End cluster for

            // Create a new cluster if this was not clustered
            if( !bHasCluster )
            {
                clusters.push_back(vector<Graph*>(1, &m_Graphs[j]));
                kCommon.push_back(m_Graphs[j]);
                if( ClusterMax.size() != clusters.size() )
                {
                    ClusterMax.push_back(1);
                    ClusterMin.push_back(1);
                }
            }
        }
    } // End graph for
  //  }


    return clusters;
}


// Calculate the probability of a graph belonging to the class with the training data specified
double GraphHandler::DeriveProbability(vector<Graph*>& kT, Graph& kInst)
{

    // 1. Generate model

    Graph kModel;
    int l=0; // Number of unique vertices


  //  vector< vector<int> > matchingTable;

  //  int iCol = 0; // Num of columns

    map<unsigned int, Vertex> kMap; // For easy retreival of vertex index for a specific ID
   // map<unsigned int, unsigned int> kMatch; // Maps matching table column -> type ID
  //  map<unsigned int, unsigned int> kTypeMatch; // Maps type ID -> table column

    int i=0;
    BOOST_FOREACH(Graph* pG, kT)
    {

        Iter_t vi, vi_end;
        Graph& G = *pG;

        // A new row in the matching table
        //matchingTable.push_back(vector<int>());

        for (tie(vi, vi_end) = vertices(G); vi != vi_end; ++vi)
        {
            unsigned int vName = get(vertex_name, G, *vi);
            if( kMap.find(vName) == kMap.end() )
            {
                Vertex vNew = add_vertex(kModel);
                //kModel[vNew].sName = m_IDTypes[vName];
                put(vertex_name, kModel, vNew, vName);

                kMap[vName] = vNew;

                l++;

                /*
                // A new column in the matching table
                for(int j=0;j<matchingTable.size();j++)
                    matchingTable[j].push_back(0); // New column
                kMatch[iCol] = vName;
                kTypeMatch[vName] = iCol;
                iCol++;*/
            }

            // Add occurance
            //matchingTable[i][kTypeMatch[vName]]++;
            i++;
        }

        graph_traits<GraphHandler::Graph>::edge_iterator firstE, lastE;
        // For each edge
        for (tie(firstE, lastE) = edges(G); firstE != lastE; ++firstE)
        {
            Vertex vSource = kMap[get(vertex_name, G, source(*firstE, G))];
            Vertex vTarget = kMap[get(vertex_name, G, target(*firstE, G))];

            pair<Edge, bool> kRet = edge( vSource, vTarget, kModel);

            if( !kRet.second )
            {
                // Add edge
                add_edge(vSource,vTarget, kModel);
            }
        }
    }


    double o = 0.0;

    vector<Graph*> kVec;
    //kVec.resize(kT.size());

    //copy(kT.begin(), kT.end(), kVec.begin());
    kVec.push_back(&kModel);
    kVec.push_back(&kInst);

    SaveGraph("model.dot", kModel );
    SaveGraph("instance.dot", kInst );

    // Let kInst -> m' be the projection of the specific instance to the model.
    // We want to count how many times this subgraph has occured in the training set.
    cout << "Before Gspan" << endl;
    pair<bool, Graph> EProj;// = GSpanStructure(kVec, 2, num_vertices(kInst));
    EProj.first = true;
    EProj.second = LargestCommonSubGraph(kModel, kInst);
    cout << "After Gspan" << endl;

    // EProj.second(m')  is the part of the model matched by E

    SaveGraph("debug.dot", EProj.second );

    if( !EProj.first )
        return 0.0;

    for(int i = 0; i<kT.size(); i++)
    {
        Graph gTmp = LargestCommonSubGraph(kInst, *kT[i]);
        o += (double)num_vertices(gTmp) / (double)min(num_vertices(kInst), num_vertices(*kT[i]));
    }

    return o / kT.size();

    for(int i = 0; i<kT.size(); i++)
    {

        /*
        kVec.clear();
        //kVec.resize(clusters[c].size());
        kVec.push_back(kT[i]);
        kVec.push_back(&EProj.second);
        pair<bool, Graph> kPTmp = GSpanStructure(kVec, 2, 0);
        if( kPTmp.first ) o += 1.0;*/

        bool bAdd = true;

        /*
        Iter_t vi, vi_end;
        for (tie(vi, vi_end) = vertices(EProj.second); vi != vi_end; ++vi)
        {
            unsigned int vName = get(vertex_name, EProj.second, *vi);
            Iter_t vi1, vi_end1;
            bool bFound = false;

            for (tie(vi1, vi_end1) = vertices(*kT[i]); vi1 != vi_end1; ++vi1)
            {
                if( get(vertex_name, *kT[i], *vi1) == vName )
                {
                    bFound = true;
                    break;
                }
            }
            if(!bFound )
            {
                bAdd = false;
                break;
            }
        }*/

        //Graph gTmp = LargestCommonSubGraph(*kT[i], EProj.second);
       // if( num_vertices(gTmp) == num_vertices(EProj.second) )
        //    o += 1.0;
        Iter_t vi, vi_end;
        double n = 0.0;
        // Go through all the vertices of N'
        for (tie(vi, vi_end) = vertices(EProj.second); vi != vi_end; ++vi)
        {

            Iter_t vi1, vi_end1;
            // Go though all the vertices of a given example of the model
            for (tie(vi1, vi_end1) = vertices(*kT[i]); vi1 != vi_end1; ++vi1)
            {
                // Intersection
                if( get(vertex_name, *kT[i], *vi1) == get(vertex_name, EProj.second, *vi) )
                    n += 1.0;
            }
        }

        o += (n * n) / (double)(num_vertices(EProj.second) * num_vertices(*kT[i]));
        //if(bAdd) o += 1.0;

    }

    //return o / (double)kT.size();

    // Calculate omax

    double oMax = 0.0;
    // sum from 1 to r (over all examples)
    for(int i = 1; i<=kT.size(); i++)
    {
        // l is the amount of unique vertices (i.e sum over all possible vertices)
        for(int n = 1; n<l; n++)
        {
            int NI = num_vertices(*kT[i-1]); // NI is cardinality of the node set of the specific example

            int iTo = min(n, NI);
            for(int e=max(1, n-l+NI); e<iTo; e++ )
            {
                double dTmp = 0.0;
                dTmp = (double)(e) / (double)(NI);
                dTmp *= (double)(e)/(double)(n);
                dTmp *= binomial(NI, e);
                dTmp *= binomial(l-NI, n-e);
                cout << "dTmp: " << dTmp << endl;
                oMax += dTmp;
            }
        }
    }

    cout << "o: " << o << " omax: " <<oMax;
    //return o / (double)kT.size();
    return o / oMax;
}

vector< Cluster* > GraphHandler::StructureClustering3(double dTheta)
{
    cout << "Structure clustering 3" << endl;
    vector<GSpanStruct> kRes = GSpan(m_Graphs, dTheta);
    cout << "Frequent subgraphs: " << kRes.size() << endl;

    for(int i=0; i<kRes.size(); i++)
    {
        Cluster* pNewC = new Cluster();
        m_pAllclusters.push_back(pNewC);
        pNewC->kMeta = kRes[i].G;

        for(int j=0; j < kRes[i].kGraphIDs.size(); j++)
            //clusters[clusters.size() - 1].push_back(&m_Graphs[kRes[i].kGraphIDs[j]]);
            pNewC->pGraphs.push_back(kRes[i].kGraphIDs[j]);//(&m_Graphs[kRes[i].kGraphIDs[j]]);

        kCommon.push_back(kRes[i].G);
    }

    cout << "Entering edit maps generation" << endl;
    GenerateFrequentEditMaps(kRes);

    return m_pAllclusters;
}


vector< Cluster* > GraphHandler::StructureClustering2()
{
    int iMaxStruct = 10;
    int iFreq = 0.1;
    //vector< vector<Graph*> > clusters;

    m_debug = false;

    vector<GSpanStruct> kRes = GSpan(m_Graphs, 0.05);
    for(int iMS = iMaxStruct; iMS >= 1; iMS--)
    {
        cout << "iMS: " << iMS << endl;

        // Should also check if the cluster we're about to add, is a subgraph of any
        // of the already existing clusters
        for(int i=0; i<kRes.size(); i++)
        {
            if(num_vertices(kRes[i].G) == iMS)
            {
                Cluster* pNewC = new Cluster();
                m_pAllclusters.push_back(pNewC);
                pNewC->kMeta = kRes[i].G;
                // Add a new cluster, first check if it is a child of any existing cluster
                //if( LargestCommonSubGraph() )
                bool bFound = false;
                for(int j=0; j < clusters.size(); j++)
                {

                    if( AddCluster(clusters[j], kRes[i].G, pNewC) )
                        bFound = true;

                }
                if(!bFound)
                {
                //clusters.push_back(vector<Graph*>());
                    clusters.push_back(pNewC);
                }

                kCommon.push_back(kRes[i].G);

                for(int j=0; j < kRes[i].kGraphIDs.size(); j++)
                    //clusters[clusters.size() - 1].push_back(&m_Graphs[kRes[i].kGraphIDs[j]]);
                    pNewC->pGraphs.push_back(kRes[i].kGraphIDs[j]);//(&m_Graphs[kRes[i].kGraphIDs[j]]);

            }
        }
        cout << "Cluster size: " << clusters.size() << endl;
    }


    return clusters;
}

bool GraphHandler::ParentsHasSubgraph(const Cluster& kC, const Graph& kG) const
{
    if( kC.pParent == 0 ) // Base case
        return false;

    // If the parent has the graph as a subgraph, return true
    if( IsSubgraphIsomorph(kG, kC.pParent->kMeta) )
        return true;

    // Interative function
    return ParentsHasSubgraph(*kC.pParent, kG );
}

bool GraphHandler::AddCluster(Cluster* pAddTo, const Graph& kMetaG, Cluster* pNewC)
{

   // if( IsSubgraphIsomorph(pAddTo->kMeta, kMetaG)
    //        && !ParentsHasSubgraph(*pAddTo, kRes[i].G) )
   // {
        // Add it as a child to this cluster or to one of its children
        // Add it to this cluster of we cannot find any children so that
        // kMetaG is contained in them.


        bool bFound = false;
        for(int i=0; i < pAddTo->pChildren.size(); i++)
        {


            if( IsSubgraphIsomorph(kMetaG, pAddTo->pChildren[i]->kMeta) )
            {
                //return AddCluster(pAddTo->pChildren[i], kMetaG, pNewC);
                if( AddCluster(pAddTo->pChildren[i], kMetaG, pNewC) )
                    bFound = true;
            }


        }
        if(bFound)
            return true;

        // We could not find any children to add it to
        if( IsSubgraphIsomorph(kMetaG,pAddTo->kMeta) )
        {
            // Do not allow children of themselves
            if( EditDistance(pNewC->kMeta, pAddTo->kMeta) != 0 )//pNewC != pAddTo )
            {
                pNewC->pParent = pAddTo;
                pAddTo->pChildren.push_back(pNewC);
                return true;
            }
            else
                return true;
        }

        return false;
  //  }
}

// Have a set of labels, V
// Generate a random graph with a fixed probability for each edge


// 1. Pick a random graph from our set of graphs
// 2. Pick an induced subgraph of this graph
// 3. Check the probabilities of the subgraph belonging to all clusters
// 4. It should be expected that it has higher probability of belonging to those clusters
// which it actually belongs.

int GraphHandler::GenerateTestExample(bool bReload, int iID, string* sVertexName)
{
    int iRn = iID;
    if( !bReload )
    {
        // 1. Select a graph at random
        if( iID == -1 )
            iRn = rand() % m_Graphs.size();
        m_pEGraph = &m_Graphs[iRn];
    }

    iExcludeGraph = iRn;

    m_kNotAdded.clear();
    m_kNotAddedEdges.clear();
    m_PGMap.clear();
    m_kPartialG = Graph();


    // Copy all the vertices into a vector and random shuffle them

    Iter_t Ibegin, Iend; tie(Ibegin, Iend) = vertices(*m_pEGraph);
    m_kNotAdded.resize(num_vertices(*m_pEGraph));

    copy(Ibegin, Iend, m_kNotAdded.begin());

    random_shuffle( m_kNotAdded.begin(), m_kNotAdded.end() );
    //---

    // Copy all the edges into a vector and random shuffle them

    graph_traits<Graph>::edge_iterator firstE, lastE;
    tie(firstE, lastE) = edges(*m_pEGraph);

    m_kNotAddedEdges.resize(num_edges(*m_pEGraph));
    copy(firstE, lastE, m_kNotAddedEdges.begin());
    random_shuffle( m_kNotAddedEdges.begin(), m_kNotAddedEdges.end() );
    //---

    // Add the initial vertex
    AddExampleVertex(true, sVertexName);

    // Now the generation of a test graph is complete
    return iRn;
}
// Checks if the two vertex sets have exactly iK elements in common
// Can be used to check for disjunction
bool GraphHandler::HasVertexCommon(const Graph& G1, const Graph& G2, int iK) const
{
    set<unsigned int> G1ID;
    Iter_t vi, vi_end;
    int n = 0;


    for (tie(vi, vi_end) = vertices(G1); vi != vi_end; ++vi)
        G1ID.insert(get(vertex_name, G1, *vi));

    for (tie(vi, vi_end) = vertices(G2); vi != vi_end; ++vi)
        if( G1ID.find(get(vertex_name, G2, *vi)) != G1ID.end() )
            n++;

    return (n == iK);
}


GraphHandler::Graph GraphHandler::GraphUnion( const vector<Graph>& kGraphs ) const
{
    list<const Graph*> kGraphsTmp;
    for( vector<Graph>::const_iterator it = kGraphs.begin(); it != kGraphs.end(); ++it )
        kGraphsTmp.push_back(&(*it));
    return GraphUnion(kGraphsTmp);
}

GraphHandler::Graph GraphHandler::GraphUnion( list<Graph>& kGraphs ) const
{
    list<const Graph*> kGraphsTmp;
    for( list<Graph>::iterator it = kGraphs.begin(); it != kGraphs.end(); ++it )
        kGraphsTmp.push_back(&(*it));
    return GraphUnion(kGraphsTmp);
}

// Returns the union of the list of graphs
GraphHandler::Graph GraphHandler::GraphUnion( list<const Graph*>& kGraphs ) const
{
    Graph kUnion;

    map<unsigned int, Vertex> vMap;

    for(list<const Graph*>::const_iterator it=kGraphs.begin(); it!=kGraphs.end(); ++it)
    {
        const Graph* G = *it;
        // 1. Add all those vertices which are not already in the union
        Iter_t vBegin, vEnd;
        for( tie(vBegin, vEnd) = vertices(*G); vBegin!=vEnd; ++vBegin )
        {
            int iId = get(vertex_name, *G, *vBegin);
            if( vMap.find(iId) == vMap.end() )
            {
                // Add it
                Vertex vNew = add_vertex(kUnion);
                //kUnion[vNew].sName = ((*G)[*vBegin]).sName;
                put(vertex_name, kUnion, vNew, iId);
                vMap[iId] = vNew;
            }
        }



        // 2. Add all the edges which are not already in the union
        graph_traits<GraphHandler::Graph>::edge_iterator firstE, lastE;
        // For each edge
        for (tie(firstE, lastE) = edges(**it); firstE != lastE; ++firstE)
        {
            Vertex vSource = vMap[get(vertex_name, **it, source(*firstE, **it))];
            Vertex vTarget = vMap[get(vertex_name, **it, target(*firstE, **it))];

            pair<Edge, bool> kRet = edge( vSource, vTarget, kUnion);

            if( !kRet.second )
                // Add edge
                add_edge(vSource,vTarget, kUnion);
        }

    }

    return kUnion;
}


template <typename Iterator>
bool next_combination(const Iterator first, Iterator k, const Iterator last)
{
   /* Credits: Mark Nelson http://marknelson.us */
   if ((first == last) || (first == k) || (last == k))
      return false;
   Iterator i1 = first;
   Iterator i2 = last;
   ++i1;
   if (last == i1)
      return false;
   i1 = last;
   --i1;
   i1 = k;
   --i2;
   while (first != i1)
   {
      if (*--i1 < *i2)
      {
         Iterator j = k;
         while (!(*i1 < *j)) ++j;
         std::iter_swap(i1,j);
         ++i1;
         ++j;
         i2 = k;
         std::rotate(i1,j,last);
         while (last != j)
         {
            ++j;
            ++i2;
         }
         std::rotate(k,i2,last);
         return true;
      }
   }
   std::rotate(first,k,last);
   return false;
}

bool GraphHandler::GraphEqual(const Graph& G1, const Graph& G2) const
{
    if( num_vertices(G1) != num_vertices(G2) )
        return false;
    if( num_edges(G1) != num_edges(G2) )
        return false;

    Iter_t G1begin, G1end; tie(G1begin, G1end) = vertices(G1);
    Iter_t G2begin, G2end; tie(G2begin, G2end) = vertices(G2);

    vector<Vertex> V1, V2;
    V1.resize(num_vertices(G1));
    V2.resize(num_vertices(G2));

    copy(G1begin, G1end, V1.begin());
    copy(G2begin, G2end, V2.begin());

    Sorter kSortG1(G1); Sorter kSortG2(G2);
    sort(V1.begin(), V1.end(), kSortG1);
    sort(V2.begin(), V2.end(), kSortG2);

    if( !std::equal(V1.begin(), V1.end(), V2.begin()) )
        return false;

    graph_traits<Graph>::edge_iterator G1Ebegin, G1Eend, G2Ebegin, G2Eend;
    tie(G1Ebegin, G1Eend) = edges(G1); tie(G2Ebegin, G2Eend) = edges(G2);

    for( tie(G1Ebegin, G1Eend) = edges(G1); G1Ebegin != G1Eend; ++G1Ebegin )
    {
        Vertex u = FindVertex(get(vertex_name, G1, source(*G1Ebegin, G1)), G2).second;
        Vertex v = FindVertex(get(vertex_name, G1, target(*G1Ebegin, G1)), G2).second;
        if( !edge(u, v, G2).second )
            return false;
    }

    return true;
}

pair<const GraphHandler::Graph*, double> GraphHandler::FindOptimalParent(const Graph& pChild,
                                    const vector<Graph>& kVec, int iIndex, const Graph& G) const
{
    Cluster* pBest = 0;
    double dFreq = -99;

    int i2 = -1;

    for(int i=0; i < m_pAllclusters.size(); i++)
    {
        if( GraphEqual(pChild, m_pAllclusters[i]->kMeta) )
        {
            i2 = FindClusterID(m_pAllclusters[i]);
            break;
        }
    }

    if( i2 == -1 )
        return make_pair((const Graph*)0,0.0);


    for(int i=0; i < m_pAllclusters.size(); i++)
    {
        // See if this cluster has pChild as a child
        if( IsSubgraphIsomorph(pChild, m_pAllclusters[i]->kMeta) )
        for(vector<Cluster*>::const_iterator it1 = m_pAllclusters[i]->pChildren.begin(); it1 != m_pAllclusters[i]->pChildren.end(); ++it1  )
        {
            if( GraphEqual((*it1)->kMeta, pChild) )
            {
                int i1 = FindClusterID((*it1));

                double dProb = GetClusterCorrelation(i1, i2);

                if( dProb >= dFreq && i1 != i2 )
                {
                    // Make sure that the new graph union contains some new information
                    vector<Graph> kTmpVec;
                    for(int j=0; j < kVec.size(); j++)
                        if( j != iIndex )
                            kTmpVec.push_back(kVec[j]);
                    kTmpVec.push_back(m_pAllclusters[i]->kMeta);
                    Graph kUnion = GraphUnion(kTmpVec);
                    if( !GraphEqual(G, kUnion) )
                    {
                        pBest = m_pAllclusters[i];
                        dFreq = dProb;
                    }
                }
            }
        }
    }

    if( pBest != 0 )
    {
        cout << "Found an optimal parent" << endl;
        cout << "Subgraph: " << IsSubgraphIsomorph(pChild, pBest->kMeta) << endl;
        return make_pair(&(pBest->kMeta), dFreq);
    }
    else
        return make_pair((const Graph*)0,0.0);
}

Cluster* GraphHandler::FindOptimalParent(const Cluster* pChild,const list<Cluster*>& kRelevantClusters) const
{
    Cluster* pBest = 0;
    unsigned int iFreq = 0;
    double dFreq = -99;


    //SaveGraph("pchild.dot", (pChild)->kMeta);
    //for(list<Cluster*>::const_iterator it = m_pAllclusters.begin(); it != m_pAllclusters.end(); ++it  )
    for(int i=0; i < m_pAllclusters.size(); i++)
    {
        // See if this cluster has pChild as a child
        for(vector<Cluster*>::const_iterator it1 = m_pAllclusters[i]->pChildren.begin(); it1 != m_pAllclusters[i]->pChildren.end(); ++it1  )
        {

            //if( EditDistance((*it1)->kMeta, pChild->kMeta) == 0 )
            if( GraphEqual((*it1)->kMeta, pChild->kMeta) )
            {
                /*
                if( (*it1)->pGraphs.size() >= iFreq )
                {
                    pBest = m_pAllclusters[i];
                    iFreq = (*it1)->pGraphs.size();
                }*/
                //if( dCorr[] )
                int i1 = FindClusterID((*it1));
                int i2 = FindClusterID(pChild);
                double dProb = GetClusterCorrelation(i1, i2);
                //double dProb = m_dCorr[i1][i2];
                if( dProb >= dFreq )
                {
                    pBest = m_pAllclusters[i];
                    dFreq = dProb;

                    //SaveGraph("child.dot", (*it1)->kMeta);
                    //SaveGraph("parent.dot", m_pAllclusters[i]->kMeta);

                }
            }

        }
    }

    return pBest;
}

// Based on the current partial graph, it performs the most likely
// action, for example adding a new vertex and connecting it, or adding
// an  between existing vertices.
GraphHandler::Graph GraphHandler::ExampleStepPredict(const Graph* pPartialGraph)
{
    // 1. Generate all subgraphs of the current partial graph
    // 2. Pick some assignement of the subgraphs to matching clusters
    // 3. Evaluate the assignment, pick the optimal one.


    // Get the clusters which have atleast one vertex in common with the partial graph

    if( pPartialGraph == 0 )
        pPartialGraph = &m_kPartialG;



    // Get all subgraphs
    vector<Graph> kTmpV; kTmpV.push_back(*pPartialGraph);
    vector<GSpanStruct> kRes = GSpan(kTmpV, 0.1);

    map<unsigned int, Cluster*> pOptCluster;
    map<Cluster*, Cluster*> pOptParent;

    // Used for saving all the clusters which are 5 edit distances away from the specific partial graph
    map< unsigned int, list<Cluster*> > kRClusterMap;

    list<Cluster*> kRelevantClusters;
    //for(list<Cluster*>::iterator it = m_pAllclusters.begin(); it != m_pAllclusters.end(); ++it  )
    for(int i=0; i < m_pAllclusters.size(); i++)
        //if( !HasVertexCommon((*it)->kMeta, *pPartialGraph, 0) ) // All those which do not have exactly 0 vertices in common
        if( !HasVertexCommon(m_pAllclusters[i]->kMeta, *pPartialGraph, 0) )
        {

            // Also have condition that the distance from any subgraph is atleast 5
            bool bSat = false;
            for(int j=0; j < kRes.size(); j++)
                //if( (EditDistance((*it)->kMeta, kRes[j].G)) < 5 )
                if( (EditDistance(m_pAllclusters[i]->kMeta, kRes[j].G)) < 5 )
                {
                    bSat = true; break;
                }
            //if( bSat ) kRelevantClusters.push_back((*it));
            if( bSat ) kRelevantClusters.push_back(m_pAllclusters[i]);


        }


    for(int j=0; j < kRes.size(); j++)
    {
        int iDistMin = 9999;
        for(list<Cluster*>::iterator it = kRelevantClusters.begin(); it != kRelevantClusters.end(); ++it  )
        {
            int iDist;
            if( (iDist = EditDistance((*it)->kMeta, kRes[j].G)) < iDistMin )
            {
                iDistMin = iDist;
                pOptCluster[j] = *it;
                pOptParent[*it] = FindOptimalParent(*it, kRelevantClusters);
            }
        }
    }

/*
    // Erase the one-vertex only subgraphs
    for(int i=0; i < kRes.size(); i++)
        if( num_vertices(kRes[i].G) < 2 )
            kRes.erase(kRes.begin() + i);
*/

    int iDistLowest = 9999;
    double dBestProb = 0.0; // The probability of the best theory so far
    m_kHypothesis.clear();
    //int iDistLowest = 0;
    list<const Graph*> kLowestHyp;

/*
    set< pair<int, int> > kChecked;

    // kHypothesis[i][j] is the jth combination of ith selection size
    vector< vector<Cluster*> > kHypothesis;
    kHypothesis.resize(kRes.size());*/

    cout << "Size: " << kRes.size() << endl;
    vector<int> kComb;
    kComb.resize(kRes.size());

    for(int j=0; j < kRes.size(); j++)
        kComb[j] = j;



    // Now, iterate through all possible combinations of subgraphs, first selection size
    // If all subgraphs are denoted G1, ..., Gn, we will form: { {G1}, {G2}, ..., {Gn}, {G1, G2}, {G1, G3}, ..., {G1, Gn} etc}
    for(int k=1; k <= kRes.size(); k++)
    {



        // Initial combination consists of (0, 1, ..., k)
        for(int j=0; j < kRes.size(); j++)
            kComb[j] = j;


        // Loop until we cannot find any more combination of sets of subgraphs with size k
        do
        {

            // In each iteration,  kComb which specifies which subgraphs
            // we are dealing with. Iterate through this specific combination
            // and check the optimal cluster for each subgraph if the combination is valid.

            // The combination is valid if:
            // (1) the partial graph is a subgraph ("is contained in") of the union of the subgraphs in the combination
            // (2) for the subgraphs in the combination, there exists another subgraph such that they have atleast one vertex in common

            // 1. Form the union and check if the partial graph is a subgraph of it
            list<const Graph*> kList;

            cout << "kComb: ";
            for(int j=0; j < k; j++)
            {
                cout << kComb[j] << " ";
                kList.push_back(&kRes[kComb[j]].G);
            }

            cout << endl;



            // For this specific combination of subgraphs, keep track of
            // which of the subgraphs have a vertex in common ("bridge")
            // This will be used for calculating probabilities
            vector<Graph> _graphs;
            for(int j=0; j < k; j++)
                _graphs.push_back(kRes[kComb[j]].G);

            list< pair<int, int> > kConnectionList;
            bool bCond2 = true;
            tie(bCond2, kConnectionList) = GetConnectionList(_graphs);





         //   m_debug = true;
            //Graph kUnion;
            Graph kUnion = GraphUnion(kList);
/*

            SaveGraph("partialgraph.dot", m_kPartialG);
            SaveGraph("uniongraph.dot", kUnion);*/

            //if( IsSubgraphIsomorph(m_kPartialG, kUnion) && (k == 1 || bCond2) )
            if( (k == 1 || bCond2) )
            {
                //cout << "Found candidate" << endl;
                int iDist = 0;
                int iSupport = 0;
                double dProb = 1.0;

                list<const Graph*> kTmpCV;
                vector<Cluster*> kClusterV;
                for(int j=0; j < k; j++)
                {
                    //Cluster* pC = GreedyFindBestCluster(kRes[kComb[j]].G, kRelevantClusters);
                    Cluster* pC = pOptCluster[kComb[j]];
                    //if( pC->pParent != 0 )
                    //if(  )
                    //{
                       // kTmpCV.push_back(&pC->pParent->kMeta);
                       // iDist += EditDistance(pC->kMeta, kRes[kComb[j]].G);
                    //}
                    if( pOptParent.find(pC) != pOptParent.end() && pOptParent[pC] != 0 )
                    {
                        kTmpCV.push_back(&pOptParent[pC]->kMeta);
                        kClusterV.push_back(pC);
                        iSupport += pOptParent[pC]->pGraphs.size();
                        iDist += EditDistance(pOptParent[pC]->kMeta, kRes[kComb[j]].G);


                        dProb *= GetProbabilityForHypothesis(pOptParent[pC]->kMeta, kRes[kComb[j]].G);
                    }
                    else
                    {
                        kTmpCV.push_back(&pC->kMeta);
                        kClusterV.push_back(pC);
                        iSupport += pC->pGraphs.size();
                        iDist += EditDistance(pC->kMeta, kRes[kComb[j]].G);

                        dProb *= GetProbabilityForHypothesis(pC->kMeta, kRes[kComb[j]].G);
                    }
                }


                // We need to get the "connection chain" of the clusters
                for( list< pair<int, int> >::const_iterator it = kConnectionList.begin(); it!=kConnectionList.end(); ++it )
                {
                    //dProb *= m_dCorr[ make_pair( kClusterV[it->first], kClusterV[it->second] ) ];
                    int i1 = FindClusterID(kClusterV[it->first]);
                    int i2 = FindClusterID(kClusterV[it->second]);
                    //dProb *= m_dCorr[ make_pair( i1, i2 ) ];
                    //dProb *= m_dCorr[i1][i2];
                    dProb *= GetClusterCorrelation(i1, i2);
                }

                cout << "dProb: " << dProb << endl;

                kUnion = GraphUnion(kTmpCV);
                int iDist2 = EditDistance(kUnion, *pPartialGraph);

                dProb *= GetProbabilityForHypothesis(kUnion, *pPartialGraph);
                //m_dCorr
                // dBestProb

                bool bIso = IsSubgraphIsomorph(*pPartialGraph, kUnion);
                bool bAlreadyExists = false;

                for( vector< pair<Graph, double> >::iterator it=m_kHypothesis.begin(); it != m_kHypothesis.end(); ++it )
                    if( EditDistance(it->first, kUnion) == 0 )
                        bAlreadyExists = true;
                if( bIso && !bAlreadyExists)
                    m_kHypothesis.push_back( make_pair(kUnion, dProb) );
                //if( iDist < iDistLowest && iDist != 0 && iDist2 != 0 && num_vertices(kUnion) > num_vertices(*pPartialGraph) )
                if( dProb > dBestProb && num_vertices(kUnion) > num_vertices(*pPartialGraph) && IsSubgraphIsomorph(*pPartialGraph, kUnion) && bIso )
                {
                    iDistLowest = iDist;
                    dBestProb = dProb;
                    kLowestHyp = kTmpCV;
                    cout << "New lowest hypothesis" << endl;

                    SaveGraph("hyp.dot", GraphUnion(kLowestHyp));
                }
            }


        } while( next_combination(kComb.begin(), kComb.begin() + k, kComb.end()) );

    }

/*
    list<const Graph*> listtmp;
    for(list<Graph>::iterator it = kLowestHyp.begin(); it!=kLowestHyp.end(); ++it)
    {
        listtmp.push_back(&(*it));
    }*/




    // Sort according to probability
    sort( m_kHypothesis.begin(), m_kHypothesis.end(), sort_pred<Graph>() );

    cout << "Probability: " << dBestProb << endl;
    if( kLowestHyp.size() != 0 )
    {
        Graph kG = GraphUnion(kLowestHyp);
        return kG;
    }
    else
    {
        return Graph();
    }

}

pair< bool, list< pair<int, int> > > GraphHandler::GetConnectionList(vector<Graph>& graphs)
{
    bool bCond2 = true;
    list< pair<int, int> > kConnectionList;
    for(int j=0; j < graphs.size(); j++)
    {

        bool bFound = false; // It is enough to find a single graph with common vertices
        for(int i=0; i < graphs.size(); i++)
        {
            // If it has no vertex in common, then set condition false and break
            if( i != j && !HasVertexCommon(graphs[j], graphs[i], 0) )
            {
                bFound = true;
                // We found that the subgraph pair (i,j) do in fact have a common vertex
                if( find(kConnectionList.begin(), kConnectionList.end(), make_pair(i, j)) == kConnectionList.end() &&
                        find(kConnectionList.begin(), kConnectionList.end(), make_pair(j, i)) == kConnectionList.end() )
                    kConnectionList.push_back( make_pair(j, i) );
                //break;
            }
        }

        if( !bFound )
        {
            bCond2 = false;
            //break;
        }
    }

    return make_pair(bCond2, kConnectionList);
}

GraphHandler::Graph GraphHandler::ExampleStepPredict4(const Graph* pPartialGraph)
{
    m_bCorrectPrediction = false;
    if( pPartialGraph == 0 )
        pPartialGraph = &m_kPartialG;
    Graph kPartial = *pPartialGraph;
    vector<Graph> C = DecomposeGraph(kPartial);

    m_iHypSupport.clear();
    vector<Graph> kGTmpVec;

    Graph kBest;
    double dBest = -9;
    Graph kBestComplete;
    int iCIndex = -1;


    for(int i=0; i < C.size(); i++)
    {
        //Cluster* pParent = FindOptimalParent(pC, kRelevantClusters);
        //pair<const Graph*, double> pParent = FindOptimalParent(C[i]);
        pair<Graph, double> kOTmp = GetOptimalParentOnline2(C[i], C, i, kPartial);
        kGTmpVec.push_back(kOTmp.first);
        m_iHypSupport.push_back(kOTmp.second);
        //editVec.push_back(kEditG);

        if( kOTmp.second > dBest )
        {
            //kHyp.push_back(*(pParent.first));
            dBest = kOTmp.second;
            kBest = kEditG;
            kBestComplete = kOTmp.first;
            iCIndex = i;
        }
        //else
        //    kHyp.push_back(C[i]);
    }

    vector<Graph> kHyp;
    if( num_vertices(kBest) != 0 )
         kHyp.push_back(kBestComplete);

    for(int i=0; i < C.size(); i++)
        if( i != iCIndex )
            kHyp.push_back(C[i]);


    if( IsSubgraphIsomorph(kBest, *m_pEGraph) && iCIndex != -1 )
        m_bCorrectPrediction = true;

    Graph kUnion = GraphUnion(kHyp);
    return kUnion;

}

GraphHandler::Graph GraphHandler::ExampleStepPredict2Preprocessed(const Graph* pPartialGraph)
{
    m_bCorrectPrediction = false;
    if( pPartialGraph == 0 )
        pPartialGraph = &m_kPartialG;
    Graph kPartial = *pPartialGraph;

    //vector<Graph> C = DecomposeGraph3(kPartial);
    // Get a vector of cluster IDs of the frequent subgraphs

    GetPatterns(kPartial);

    //m_iIDVector now contains the ID of the frequent subgraphs which are contained in this partial graph
    EditOperationMap kRes;
    cout << "Before sum, size: " <<  m_iIDVector.size() <<  endl;
    for(int i=0; i < m_iIDVector.size(); i++)
        kRes += m_kEditMaps[m_iIDVector[i]];
    cout << "After sum" << endl;

    if( m_iIDVector.size() == 0 )
    {
        // If we cannot find any frequent subgraphs
        kRes = EditOperationMap(*pPartialGraph);
    }

    cout << "Before get optimal" << endl;
    // Now get the valid edit operations
    pair< EditOperationMap::VEditOp, int > kOptV = kRes.GetOptimalVertexEditOp(kPartial);
    pair< EditOperationMap::VEditOp, int > kOptE = kRes.GetOptimalEdgeEditOp(kPartial);
    cout << "After get optimal" << endl;

    pair<Graph, Graph> kResultGraphs;
    int iBestSupport = -1;

    if( kOptV.second > kOptE.second )
    {
        kResultGraphs = EditOperationMap::FormResultOfVertexEditOp(kPartial, kOptV.first);
        iBestSupport = kOptV.second;
    }
    else
    {
        kResultGraphs = EditOperationMap::FormResultOfEdgeEditOp(kPartial, kOptE.first);
        iBestSupport = kOptE.second;
    }

    if( IsSubgraphIsomorph(kResultGraphs.second, *m_pEGraph) && iBestSupport != -1 )
        m_bCorrectPrediction = true;

    return kResultGraphs.first;
}

GraphHandler::Graph GraphHandler::ExampleStepPredict2(const Graph* pPartialGraph)
{
   // cout << "Entering predict" << endl;
    //if( pPartialGraph != 0 )
    //SaveGraph("predictg.dot", *pPartialGraph);

    m_bCorrectPrediction = false;
    if( pPartialGraph == 0 )
        pPartialGraph = &m_kPartialG;

    // 1. Generate relevant clusters (all for now)
    list<Cluster*> kRelevantClusters;
    for(int i=0; i < m_pAllclusters.size(); i++)
        kRelevantClusters.push_back(m_pAllclusters[i]);
    // -------------------------


    Graph kPartial = *pPartialGraph;
    //vector<Graph> C = MetisPartition(kPartial);//RareAlgorithm(kPartial);
    boost::timer timer;

    //cout << "Decompose graph" << endl;
    vector<Graph> C = DecomposeGraph3(kPartial);
    //vector<Graph> C = DecomposeGraph3(kPartial);
   // cout << "Finished Decompose graph" << endl;
    //vector<Graph> C = MetisPartition(kPartial);
    //vector<Graph> C; C.push_back(*pPartialGraph);
    //std::cout << "Decompose graph: " <<  timer.elapsed() << " seconds." << std::endl;

    m_kPSubgraphs = C;

    list< pair<int, int> > kConnectionList;
    bool bCond = true;

    //tie(bCond, kConnectionList) = GetConnectionList(C);

    list<Graph> kHyp;


    //cout << "C size: " << C.size() << endl;
    /*
    for(int i=0; i < C.size(); i++)
    {
        Cluster* pC = GreedyFindBestCluster(C[i], kRelevantClusters);

        if( pC != 0 )
        {
            Cluster* pParent = FindOptimalParent(pC, kRelevantClusters);
            if( pParent != 0 )
                kHyp.push_back(pParent->kMeta);
            else
                kHyp.push_back(pC->kMeta);
        }
    }*/

    timer.restart();
    double dBest = -9.0;
    //const Graph* pBest = 0;
    Graph pBest;
    int iCIndex = -1;
    //vector<Graph> editVec;
    vector< pair<Graph,int> > editVec;

/*
FindOptimalParent(const Graph& pChild,
                                        const vector<Graph>& kVec, int iIndex, const Graph& G
  */

    m_iHypSupport.clear();
    vector<Graph> kGTmpVec;
    /*
    for(int i=0; i < C.size(); i++)
    {
        //Cluster* pParent = FindOptimalParent(pC, kRelevantClusters);
        //pair<const Graph*, double> pParent = FindOptimalParent(C[i]);
        pair<Graph, int> kOTmp = GetOptimalParentOnline(C[i], C, i, kPartial);
        kGTmpVec.push_back(kOTmp.first);
        m_iHypSupport.push_back(kOTmp.second);
        editVec.push_back(kEditG);


        pair<const Graph*, double> pParent;// = FindOptimalParent(C[i], C, i, kPartial);
        pParent.first = &kOTmp.first; pParent.second = (double)kOTmp.second;
        if( pParent.first != 0 && pParent.second > dBest )
        {
            //kHyp.push_back(*(pParent.first));
            dBest = pParent.second;
            pBest = *pParent.first;
            iCIndex = i;
        }
        //else
        //    kHyp.push_back(C[i]);
    }*/

    pair<Graph, int> kBest;
    kBest.second = -1;



    //cout << "C size:" << C.size() << endl;
    for(int i=0; i < C.size(); i++)
    {
        //SaveGraph("tmp1.dot", C[i]);
       // cout << "Before get opti parent " << i << endl;
        pair<Graph, int> kOTmp = GetOptimalParentOnline(C[i], C, i, kPartial);
      //  cout << "After get opti parent " << i << endl;
        m_iHypSupport.push_back(kOTmp.second);
        kGTmpVec.push_back(kOTmp.first);
        bool bFound = false;
        int j=0;
        /*
        for(j=0; j < editVec.size(); j++)
        {
            //SaveGraph("tmp.dot", kEditG);
            //SaveGraph("tmp2.dot", editVec[j].first);
            //cout << "GE: " << GraphEqual(kEditG, editVec[j].first) << endl;
            //cout << "Edit distance: " << GraphEqual(kEditG, editVec[j].first) << endl;
            if( GraphEqual(kEditG, editVec[j].first) )
            {
                //cout << "Equal" << endl;
                bFound = true;
                editVec[j].second += kOTmp.second;
                break;
            }
        }*/

        if( !bFound )
        {editVec.push_back( make_pair(kEditG, kOTmp.second) ); j=editVec.size() - 1; }

        if( editVec[j].second > kBest.second )
        {
            kBest = editVec[j];
            iCIndex = i;
        }
    }


    m_kPSubgraphs = kGTmpVec;


    //if( num_vertices(pBest) != 0 )
     //   kHyp.push_back(pBest);

    if( num_vertices(kBest.first) != 0 )
         //kHyp.push_back(kBest.first);
        kHyp.push_back(kGTmpVec[iCIndex]);
    for(int i=0; i < C.size(); i++)
        if( i != iCIndex )
            kHyp.push_back(C[i]);

    //SaveGraph("partial.dot", kPartial);
    //SaveGraph("tmp.dot", editVec[iCIndex]);

    /*
    cout << "iCIndex: " << iCIndex << endl;
    cout << "c total: " << C.size() << endl;

    for(int i=0; i < editVec.size(); i++)
    {
        stringstream sFilename;
        sFilename << "./debuggraphs/editvec" << i << ".dot";
        SaveGraph(sFilename.str(), editVec[i]);
    }*/

    /*
    cout << "iCIndex: " << iCIndex << endl;
    cout << "c total: " << C.size() << endl;
    cout << "editVec: " << editVec.size() << endl;*/
    if( IsSubgraphIsomorph(kBest.first, *m_pEGraph) && kBest.second != -1 )
        m_bCorrectPrediction = true;




    //std::cout << "Find optimal parent: " <<  timer.elapsed() << " seconds." << std::endl;


    Graph kUnion = GraphUnion(kHyp);
    return kUnion;

}

vector<GraphHandler::Graph> GraphHandler::MetisPartition(const Graph& G)
{
    SaveMetisGraph("tmpmetis", G);

    int iParts = (int) ceil( (double)num_vertices(G) / 3.0 );

    if( iParts <= 1 )
    {
        vector<Graph> kTmp;
        kTmp.push_back(G);
        return kTmp;
    }

    stringstream sTmp;
    sTmp << "../thirdparty/metis-4.0.3/kmetis tmpmetis" << " " <<iParts;

    system (sTmp.str().c_str());


    stringstream sResult;
    sResult << "tmpmetis.part." <<  iParts;
    return LoadMetisPartition(sResult.str().c_str(), G);
}

GraphHandler::Graph GraphHandler::ExampleStepPredict3(const Graph* pPartialGraph)
{
    // 1. For each partition of the vertex set, form the induced subgraphs
    // 2. For each induced subgraph, find optimal parent
    // 3. Form union of these, pick the best hypothesis

    if( pPartialGraph == 0 )
        pPartialGraph = &m_kPartialG;
    Graph kPartial = *pPartialGraph;



    /*
    vector<int> ver;
    //Iter_t Ibegin, Iend; tie(Ibegin, Iend) = vertices(kPartial);
    ver.resize(num_vertices(kPartial));
    for(int i=0; i < num_vertices(kPartial); i++)
        ver[i] = i;


    //copy(Ibegin, Iend, ver.begin());

    try {
        SetPartition::iterator_k it(num_vertices(kPartial), min(3, (int)num_vertices(kPartial)));

      while (true) {
        std::cout << it << " : " << it.subsets() << " : ";

        std::auto_ptr<std::vector<std::vector<int> > >
          part = it[ver];
        std::cout << *part << '\n';
        cout.flush();

        ++it;
      }
    } catch (std::overflow_error&) { return Graph(); }
*/
    return Graph();
}

// For a graph, finds the cluster which is the closest to it in terms
// of the graph edit distance, if there are several, it picks the one
// with is more frequently occuring in the training set (should thus be more probable)
Cluster* GraphHandler::GreedyFindBestCluster(const Graph& G, const list<Cluster*>& kTmpV)
{
    Cluster* pBestC = 0;
    double dBest = 9999;

    int dFreq = 0;

    //for(list<Cluster*>::const_iterator it = kTmpV.begin(); it != kTmpV.end(); ++it  )
    for( vector<Cluster*>::const_iterator it = m_pAllclusters.begin(); it != m_pAllclusters.end(); ++it )
    {
        // More probable with the closer clusters
        int iDist = EditDistance(G, (*it)->kMeta);

        // The proportion of the cluster of all graphs * 1 / the distance
        double dValue;
        /*
        if( iDist != 0)
            dValue = 1.0 / (double)iDist * ((double)clusters[i]->pGraphs.size() / (double)m_Graphs.size());
        else
            dValue = (double)clusters[i]->pGraphs.size() / (double)m_Graphs.size();*/

        dValue = (double)iDist;

        if( dValue <= dBest )//&& IsSubgraphIsomorph(G, (*it)->kMeta) )
        {
            if( dValue < dBest )
            {
                dBest = dValue;
                dFreq = (*it)->pGraphs.size();
                pBestC = (*it);
            }
            else
            {
                // Equal, choose the most frequent
                if( (*it)->pGraphs.size() > dFreq )
                {
                    dBest = dValue;
                    dFreq = (*it)->pGraphs.size();
                    pBestC = (*it);
                }
            }
        }
    }
/*
    if( pBestC->pParent != 0 )
    {
        return pBestC->pParent->kMeta;
    }
    else
        return pBestC->kMeta;*/
    return pBestC;
}

void GraphHandler::ExampleStepNext()
{
    // Add a new vertex and connect it or add an edge that's not already occuring
    //iRn = rand() % num_vertices(*m_pEGraph);

    // Count the number of available vertex insertions
    int iNum = 0;
    /*
    for (tie(vi, vi_end) = vertices(*m_pEGraph); vi != vi_end; ++vi, i++)
        if( !FindVertex(get(vertex_name, *m_pEGraph, *vi), m_kPartialG).first )
            iNum++;*/


    int iVertexC = VertexAdditionsPossible();
    iNum += iVertexC;

    for(int i=0; i < m_kNotAddedEdges.size(); i++)
    {
        Edge E = m_kNotAddedEdges[i];
        Vertex vS = source(E, *m_pEGraph);
        Vertex vT = target(E, *m_pEGraph);

        // If both vertices of the edge exist, it is a possible edge addition
        if( m_PGMap.find(vS) != m_PGMap.end() && m_PGMap.find(vT) != m_PGMap.end() )
            iNum++;
    }

    if(iNum == 0)
        return;


    int iRn = rand() % iNum;
    if( iRn < iVertexC )
        // Add a new vertex from m_pEGraph and an edge to it
        // The new vertex must be connected to some of the already existing vertices
        AddExampleVertex(false);
    else
        // Add a random edge
        // A valid random edge is one such that both vertices are existing as well as
        // the edge itself exists in the graph
        AddExampleEdge();


}

// Given the current partial graph, how many vertex additions are possible?
int GraphHandler::VertexAdditionsPossible()
{
    int n = 0;
    for(int i=0; i < m_kNotAdded.size(); i++)
    {
        bool bValid = false;
        Vertex VAdd = m_kNotAdded[i];

        graph_traits<Graph>::out_edge_iterator firstE, lastE;
        for (tie(firstE, lastE) = out_edges(VAdd, *m_pEGraph); firstE != lastE; ++firstE)
            if( source(*firstE,*m_pEGraph) == VAdd )
                // So the target must be a vertex already existing
                if( m_PGMap.find(target(*firstE,*m_pEGraph)) != m_PGMap.end() )
                {
                    bValid = true;
                    break;
                }

        if( bValid )
            n++;
    }

    return n;
}

void GraphHandler::AddExampleVertex(bool bForceAdd, string* sVertexName)
{

    // Retreive a vertex not yet added, and which is possible to add at this stage


    for(int i=0; i < m_kNotAdded.size(); i++)
    {
        // Check if this vertex is addable at this time.
        // This means that there is some edge in the example graph between some already added vertex
        // in the partial graph (so we don't get disconnected parts).

        Vertex VAdd;
        if( sVertexName == 0 )
            VAdd = m_kNotAdded[i];
        else
        {
            if( m_TypesID.find(*sVertexName) == m_TypesID.end() )
                return;
            for(int j=0; j < m_kNotAdded.size(); j++)
                if( get(vertex_name, *m_pEGraph, m_kNotAdded[j]) == m_TypesID[*sVertexName] )
                    VAdd = m_kNotAdded[j];
        }
        bool bValid = false;

        graph_traits<Graph>::out_edge_iterator firstE, lastE;
        Vertex vAlreadyExistV;
        Edge kToRemove;

        for (tie(firstE, lastE) = out_edges(VAdd, *m_pEGraph); firstE != lastE; ++firstE)
        {
            if( source(*firstE,*m_pEGraph) == VAdd )
            {
                // So the target must be a vertex already existing
                if( m_PGMap.find(target(*firstE,*m_pEGraph)) != m_PGMap.end() )
                {
                    bValid = true;
                    vAlreadyExistV = m_PGMap[target(*firstE,*m_pEGraph)];
                    kToRemove = *firstE;
                    break;
                }
            }
        }

        // Force add is used to add the first vertex for instance
        if( bForceAdd || bValid )
        {

            // Add a new vertex to the partial graph and copy all the data from the example graph
            Vertex V = add_vertex(m_kPartialG);
            //m_kPartialG[V].sName = (*m_pEGraph)[VAdd].sName;
            put(vertex_name, m_kPartialG, V, get(vertex_name, *m_pEGraph, VAdd));

            // Connect it as well using the edge we found
            if( bValid )
            {
                add_edge(V, vAlreadyExistV, m_kPartialG);
                m_kNotAddedEdges.erase(find( m_kNotAddedEdges.begin(), m_kNotAddedEdges.end(), kToRemove));

                //cout << "Adding new vertex " << m_kPartialG[V].sName << " connected from " << m_kPartialG[vAlreadyExistV].sName << endl;
            }

            // Save the correspondence vertex of m_pEGraph -> vertex of partial graph
            m_PGMap[VAdd] = V;

            // Remove it from m_kNotAdded
            m_kNotAdded.erase(m_kNotAdded.begin()+i);

           return;
        }
    }


}

void GraphHandler::AddExampleEdge()
{

    for(int i=0; i < m_kNotAddedEdges.size(); i++)
    {
        Edge E = m_kNotAddedEdges[i];
        Vertex vS = source(E, *m_pEGraph);
        Vertex vT = target(E, *m_pEGraph);

        // If both vertices of the edge exist
        if( m_PGMap.find(vS) != m_PGMap.end() && m_PGMap.find(vT) != m_PGMap.end() )
        {
            add_edge(m_PGMap[vS], m_PGMap[vT], m_kPartialG);
            //cout << "Adding edge between existing vertices: " << m_kPartialG[m_PGMap[vS]].sName <<
           //         ", " << m_kPartialG[m_PGMap[vT]].sName << endl;
            m_kNotAddedEdges.erase(m_kNotAddedEdges.begin()+i);

            return;
        }
    }
}

double GraphHandler::GetClusterCorrelation(int i, int j) const
{
    double dSupG1 = (double)(m_pAllclusters[i])->pGraphs.size();
    double dSupG2 = (double)(m_pAllclusters[j])->pGraphs.size();

    double dSupG1G2 = JointGraphSupp(m_pAllclusters[i], m_pAllclusters[j]);// / (double)m_Graphs.size();

    //return dSupG1G2 / (dSupG1+dSupG2);
    return (dSupG1G2 - dSupG1*dSupG2) / sqrt(dSupG1*dSupG2*(1-dSupG1)*(1-dSupG2));;
}

void GraphHandler::CalculateClusterCorrelations()
{
    m_dCorr.resize(m_pAllclusters.size());
    // For each pair of cluster (Ci, Cj), calculate their correlation
    //for(list<Cluster*>::iterator it = m_pAllclusters.begin(); it != m_pAllclusters.end(); ++it  )
    for(int i=0; i < m_pAllclusters.size(); i++)
    {
        m_dCorr[i].resize(m_pAllclusters.size());
        //cout << j << "/" << m_pAllclusters.size() << " (" << (double)j / (double)m_pAllclusters.size() << "%\n";
        cout << i << "/" << m_pAllclusters.size() << " (" << (double)i / (double)m_pAllclusters.size() << "%\n";
        for(int j=0; j < m_pAllclusters.size(); j++)
        {
            //if( *it != *it1 )
            if( i != j )
            {
                //CPair kPair = make_pair(*it, *it1);
                CPair kPair = make_pair(i, j);
                //bool b1 = m_dCorr.find(kPair) == m_dCorr.end();
                //bool b2 = m_dCorr.find(make_pair(*it1, *it)) == m_dCorr.end();
                //bool b2 = m_dCorr.find(make_pair(i, j)) == m_dCorr.end();
                // If this combination has not yet been recorded
                bool b1 = true;
                bool b2 = true;
                if( b1 && b2 )
                {
                    // Calculate Pearson’s Correlation Coefficient, from article "Efficient Discovery of Frequent Correlated Subgraph Pairs"
                    // theta(g1, g2)= ( supp(g1,g2) - supp(g1)*supp(g2) ) / sqrt( supp(g1)*supp(g2)(1-supp(g1))(1-supp(g2)) )
                    //double dSupG1 = (double)(*it)->pGraphs.size(); /// (double)m_Graphs.size();
                    //double dSupG2 = (double)(*it)->pGraphs.size(); /// (double)m_Graphs.size();
                    double dSupG1 = (double)(m_pAllclusters[i])->pGraphs.size();
                    double dSupG2 = (double)(m_pAllclusters[j])->pGraphs.size();
                    //double dSupG1G2 = JointGraphSupp(*it, *it1);// / (double)m_Graphs.size();
                    double dSupG1G2 = JointGraphSupp(m_pAllclusters[i], m_pAllclusters[j]);// / (double)m_Graphs.size();
                    //double dTheta = (dSupG1G2 - dSupG1*dSupG2) / sqrt(dSupG1*dSupG2*(1-dSupG1)*(1-dSupG2));
                    //m_dCorr[kPair] = dTheta;
                    // The number of times we've observed both of them divided by the number of times either of them have been in a sample
                    //m_dCorr[kPair] = dSupG1G2 / (dSupG1+dSupG2);
                    m_dCorr[i][j] = dSupG1G2 / (dSupG1+dSupG2);;
                    //m_dCorr[make_pair(*it1, *it)] = m_dCorr[kPair];
                    //m_dCorr[make_pair(i, j)] = m_dCorr[kPair];
                    //m_dCorr[j][i] = m_dCorr[i][j];
                }
            }
            else
                //m_dCorr[make_pair(i, j)] = 1.0;
                m_dCorr[i][j] = 1.0;
        }
    }
}

double GraphHandler::JointGraphSupp(const Cluster* C1, const Cluster* C2) const
{
    //static vector<Graph*> kTmp;
    static vector<int> kTmp;
    if( kTmp.capacity() == 0 )
    {
        kTmp.reserve(5000);
    }
    kTmp.clear();

    set_intersection(C1->pGraphs.begin(),C1->pGraphs.end(), C2->pGraphs.begin(), C2->pGraphs.end(), back_inserter(kTmp));

    return (double)kTmp.size();
}



double GraphHandler::GetProbabilityForHypothesis(const Graph& gHyp, const Graph& gPartial)
{
    // 1. Determine which vertices (and edges) in the partial graph
    //    are connected to the vertices which are only in the hypothesis (and not in the partial)
    // 2. Create the graph which consists of those vertices as well as the H\PG graph



    //SaveGraph("diffg.dot", gHyp, &kDiffG);
    //SaveGraph("hypi.dot", gHyp);
    //SaveGraph("diffg.dot", kDiffG);
    //SaveGraph("partial.dot", gPartial);

    // With the "difference graph" in place, calculate how many times it has been observed divided by
    // the total number of examples

    Graph kDiffG = GetDifferenceGraph(gHyp, gPartial);

    int j = 0;
    for(int i=0; i < m_Graphs.size(); i++)
    {
        if( IsSubgraphIsomorph(kDiffG, m_Graphs[i]) )
            j++;
    }

    return (double)j / (double)m_Graphs.size();
}

// This functions does the following:
// 1. It picks those vertices that are in the hypothesis and not in the partial
// That is, gHyp \ gPartial
// 2. It then also adds the vertices from gPartial which are connected to some vertex in gHyp
GraphHandler::Graph GraphHandler::GetDifferenceGraph(const Graph& gHyp, const Graph& gPartial)
{
    Iter_t vi, vi_end, next;
    Graph kDiffG;
    map<Vertex, Vertex> kMap;

    // We must first create the graph of the hypothesis with the partial graph removed from it
    for (tie(vi, vi_end) = vertices(gHyp); vi != vi_end; ++vi)
    {
        // If the vertex does not exist in the partial graph, add it to kDiffG
        if( !FindVertex(get(vertex_name, gHyp, *vi), gPartial).first )
        {
            // put(property, g, x, v)
            Vertex vNew = add_vertex(kDiffG);
            put(vertex_name, kDiffG, vNew,get(vertex_name, gHyp, *vi));
            kMap[*vi] = vNew;
        }
    }

    list<Vertex> toAdd;
    list< pair<Vertex, Vertex> > toAddEdge;

    // Add the internal edges
    // Examine edges and add those which the source and target is in the kDiffG
    graph_traits<Graph>::edge_iterator firstE, lastE;
    for (tie(firstE, lastE) = edges(gHyp); firstE != lastE; ++firstE)
    {
        bool bSFound = kMap.find(source(*firstE, gHyp)) != kMap.end();
        bool bTFound = kMap.find(target(*firstE, gHyp)) != kMap.end();
        if( bSFound && bTFound )
        {
            // Add the edge, both are in the difference graph
            add_edge(kMap[source(*firstE, gHyp)], kMap[target(*firstE, gHyp)], kDiffG);
        }
        else if( bSFound && !bTFound )
        {
            // So only one of the vertices are in the difference graph, target is not
            // We should add the target vertex later as well as the edge, record it
            if( find(toAdd.begin(), toAdd.end(), target(*firstE, gHyp)) == toAdd.end() )
            {
                toAdd.push_back(target(*firstE, gHyp));
            }
            toAddEdge.push_back( make_pair(target(*firstE, gHyp), kMap[source(*firstE, gHyp)]) );
        }
        else if( !bSFound && bTFound )
        {
            // Same as above
            if( find(toAdd.begin(), toAdd.end(), source(*firstE, gHyp)) == toAdd.end() )
            {
                toAdd.push_back(source(*firstE, gHyp));
            }
            toAddEdge.push_back( make_pair(source(*firstE, gHyp), kMap[target(*firstE, gHyp)]) );
        }
    }

    for( list<Vertex>::iterator it=toAdd.begin(); it != toAdd.end(); ++it )
    {
        Vertex vNew = add_vertex(kDiffG);
        put(vertex_name, kDiffG, vNew,get(vertex_name, gHyp, *it));
        kMap[*it] = vNew;
    }

    for( list< pair<Vertex, Vertex> >::iterator it=toAddEdge.begin(); it != toAddEdge.end(); ++it )
    {
        add_edge(kMap[it->first], it->second, kDiffG);
    }

    return kDiffG;
}

int GraphHandler::GetGraphID(const Graph* pG) const
{
    for(int i=0; i < m_Graphs.size(); i++)
        if(&m_Graphs[i] == pG)
            return i;
    return -1;
}

template<class Archive>
void GraphHandler::serialize(Archive &ar, const unsigned int version)
{
    ar & m_Graphs;
    /*
    ar & m_dCorr;
    ar & clusters;
    ar & m_pAllclusters;*/

    //ar & m_dCorr;
    ar & clusters;
    ar & m_pAllclusters;

    ar & m_IDTypes;
    ar & m_TypesID;
    ar & m_Types;
}

void GraphHandler::Save(string sFilename)
{
    ofstream ofs(sFilename.c_str(), ios::binary);

    boost::archive::binary_oarchive oa(ofs);
    oa << (*this);
    ofs.close();
}

void GraphHandler::Load(string sFilename)
{
    std::ifstream ifs(sFilename.c_str(), ios::binary);


    boost::archive::binary_iarchive oa(ifs);
    oa >> (*this);

    ifs.close();
}

int GraphHandler::FindClusterID(const Cluster* pC) const
{
    for(int i=0; i < m_pAllclusters.size(); i++)
    {
        if( m_pAllclusters[i] == pC )
            return i;
    }
    return -1;
}

void PrintVector(vector<GraphHandler::Vertex> v)
{
    for(int i=0; i < v.size(); i++)
    {
        cout << v[i] << " ";
    }
    cout << endl;
}

bool GraphHandler::CheckPartition(const vector<int>& kPart, int iG) const
{
    for(int i=0; i < kPart.size(); i++)
        if( kPart[i] < 2 || kPart[i] >= iG )
            return false;
    return true;
}

vector<GraphHandler::Graph> GraphHandler::ConnectedSubgraphs(const Graph& G) const
{
    int iSeparations = 1;
    Partition kPartition(num_vertices(G) + iSeparations); // One vertex in common

    // Minimally they shall share 1 common vertex
    // Maximally each pair has a common vertex

    vector< vector< vector<Vertex> > > kRes;


    cout << "num vertices: " << num_vertices(G) << endl;

    vector<Vertex> vVec;
    Iter_t Ibegin, Iend; tie(Ibegin, Iend) = vertices(G);
    vVec.resize(num_vertices(G));

    // The vector of vertices
    copy(Ibegin, Iend, vVec.begin());



    for( Partition::Iterator it = kPartition.begin(); it != kPartition.end(); ++it )
    {
        // For a specific partition of the number of vertices, i.e |G1|+...+|Gk|=|G|

        cout << "partition: " ;
        for(int i=0; i < (*it).size(); i++)
        {
            cout << (*it)[i] << " ";
        }

        //if( (*it).size() == 1 )
         //   (*it)[] -= iSeparations; // No separations here

        // If we are demanding x separations, we must demand to have atleast x(+1?) partitions

        cout << endl;
        cout << endl;

        if( iSeparations < (*it).size() )
        {

        vector<Vertex> kComb;
        kComb.resize(num_vertices(G));
        copy(Ibegin, Iend, kComb.begin());


        // (*it)[i] is  |Gi|
        // This generates all possible combinations of selecting |G1| vertices from our vertex set
        do
        {
            if( CheckPartition(*it, num_vertices(G)) )
            {
            // For a specific combination of G1 vertices selected, we want to generate G2
            // It must have atleast 1 vertex in common with some Gj, j < i.

            cout << "kComb: ";
            for(int j=0; j < (*it)[0]; j++)
            {
                cout << kComb[j] << " ";
            }

            cout << endl;

            vector<Vertex>::iterator kCEnd = kComb.begin() + (*it)[0];

            vector<Vertex> vTmp;
            vTmp.reserve(max(vVec.size(), kComb.size()));
            //The difference of two sets is formed by the elements that are present in the first set, but not in the second one.
            set_difference (vVec.begin(), vVec.end(), kComb.begin(), kCEnd, back_inserter( vTmp ));

            // All of the added vertices may be used as common vertices
            vector<Vertex> vCTmp;
            vCTmp.resize((*it)[0]);
            copy(kComb.begin(), kCEnd, vCTmp.begin());


            vector< vector< vector<Vertex> > > kResTmp;
            bool b1;
            cout << "Outer call" << endl;
            tie(kResTmp, b1) = ConnectedSubgraphInternal(vTmp, vCTmp, 1, *it, iSeparations, 0);
            if( b1 )
            {
                cout << "Returned true" << endl;
                for( int k = 0; k < kResTmp.size(); k++ )
                {
                    // For each of these subsequent sets of vertices, add it to the result after combining it
                    // with kComb
                    //kRes
                    vector< vector<Vertex> > kNewComb = kResTmp[k];

                    // kTotal defines the new graph we're adding
                    //vector< Vertex > kTotal = kComb;
                   // kTotal.push_back(vCommon[i]);



                    // This new combination is this new graph + all possible combinations of subsequent graphs retreived by
                    // the recursive call
                    vector<Vertex> kCombTmp;
                    kCombTmp.resize((*it)[0]);
                    copy(kComb.begin(), kComb.begin() + (*it)[0], kCombTmp.begin());
                    kNewComb.push_back( kCombTmp );

                    cout << "----------FINAL SOLUTION----------:" << endl;
                    cout << "Partitions: " << kNewComb.size() << endl;
                    for(int j=0; j < kNewComb.size(); j++)
                    {
                        cout << "{ ";
                        PrintVector(kNewComb[j]);
                        cout << " }, \n";
                    }

                    kRes.push_back( kNewComb );

                }
            }
            else
            {
                cout << "Returned false" << endl;
                vector< vector<Vertex> > kNewComb;
                kNewComb.push_back( kComb );
                kRes.push_back( kNewComb );
            }
            }
        } while( next_combination(kComb.begin(), kComb.begin() + (*it)[0], kComb.end()) );

        }
    }

    cout << "--------------------------" << endl;
    return vector<GraphHandler::Graph>();
}

// This generates
// iK is the number of vertices which this set should contain
// The common vector are those vertices which have already been added
// Any of these may be picked for this new set of vertices as the common vertex



// Returns: A vector of possible combinations, where each combination is a vector where each element is a vector of vertices
// for that particular graph
pair< vector< vector< vector<GraphHandler::Vertex> > >, bool>
    GraphHandler::ConnectedSubgraphInternal(vector<Vertex> vVec, vector<Vertex> vCommon, int iN, vector<int>& kPartition, int iNSeparatorsTotal, int iSep) const
{
    vector< vector< vector<Vertex> > > kRes;

    cout << "vVec: ";
    PrintVector(vVec);

    cout << "vCommon: ";
    PrintVector(vCommon);



    cout << "iN: " << iN << ", kPart[" << iN << "]=" << kPartition[iN] << endl;


    if( iN >= kPartition.size() )
        return make_pair( vector< vector< vector<Vertex> > >(), false);

    int iK = kPartition[iN];
    // No more vertices to add
    if( vVec.size() == 0 )
        return make_pair( vector< vector< vector<Vertex> > >(), false);



    // If there are not enough vertices left, return false
    if( 1+vVec.size() < iK )
        return make_pair( vector< vector< vector<Vertex> > >(), false);

    if( iK == 1 )
        return make_pair( vector< vector< vector<Vertex> > >(), false);

    // 1. Let m=|vCommon|, there are m different possible choices of a common vertex with the previous
    // sets of vertices. For each of these
    for(int i=0; i < vCommon.size(); i++)
    {
        // We can now select all possible combinations of vertices from those left in vVec (none of those in
        // common are in vVec so these are fully unused vertices).

        vector<Vertex> kComb;
        kComb.resize(vVec.size());
        copy(vVec.begin(), vVec.end(), kComb.begin());


        cout << "vCommon[" << i << "]=" << vCommon[i] << endl;

        int iC = 0;
        //do
        //{
            int iSepTmp = iSep;
            // For a specific combination, get the subsequent sets by recursive calls

            vector<Vertex> kCombTmp;
            kCombTmp.resize((iK-1));
            copy(kComb.begin(), kComb.begin() + (iK-1), kCombTmp.begin());

            // We want to remove those vertices which are in kComb from vVec and also
            // add the selected common vertex to vCommon
            vector<Vertex>::iterator kCEnd = kComb.begin() + (iK-1);

            vector<Vertex> vTmp;
            vTmp.reserve(max(vVec.size(), kComb.size()));
            //The difference of two sets is formed by the elements that are present in the first set, but not in the second one.
            set_difference (vVec.begin(), vVec.end(), kComb.begin(), kCEnd, back_inserter(vTmp));
            //vTmp.push_back(vCommon[i]);

            //assert( vTmp.size() == iK );

            // Form a new common vector, add all the old common and this particular combination
            vector<Vertex> vCTmp;
            vCTmp.resize(vCommon.size());
            copy(vCommon.begin(), vCommon.end(), vCTmp.begin()); // Copy from vCommon->vCTmp

            cout << "vTmp: ";
            PrintVector(vTmp);

            cout << "vCTmp: ";
            PrintVector(vCTmp);

            // if iSep == iNSeparatorsTotal, then we are not allowed to add any more separator candidates
            // We must pick one of the already added (we must always do this, this is done by the outer loop)
            // Otherwise we must add a new separator (as well as pick an old one)

            if( iSep != iNSeparatorsTotal )
            {
                vCTmp.resize(vCommon.size() + (iK-1));
                cout << "Adding kComb to common list" << endl;
                copy(kComb.begin(), kCEnd, vCTmp.begin() + vCommon.size()); // Copy from kComb->vCTmp
                iSepTmp++;
                cout << "vCTmp: ";
                PrintVector(vCTmp);
            }

            vector< vector< vector<Vertex> > > kResTmp;
            bool b1 = true;

            if( (iN + 1) != kPartition.size() )
            {
                //iN++;
                cout << "Inner call" << endl;
                tie(kResTmp, b1) = ConnectedSubgraphInternal(vTmp, vCTmp, iN + 1, kPartition, iNSeparatorsTotal, iSepTmp );
            }
            else
            {
                // kTotal defines the new graph we're adding
                vector< vector<Vertex> > kNewComb;
                vector< Vertex > kTotal = kCombTmp;
                kTotal.push_back(vCommon[i]);
                kNewComb.push_back(kTotal);
                kRes.push_back( kNewComb );

                cout << "Adding :";
                PrintVector(kTotal);
            }

            if( !b1 )
                 return make_pair( vector< vector< vector<Vertex> > >(), false);

            // For each combination
            for( int k = 0; k < kResTmp.size(); k++ )
            {

                vector< vector<Vertex> > kNewComb = kResTmp[k];


                // kTotal defines the new graph we're adding
                vector< Vertex > kTotal = kCombTmp;
                kTotal.push_back(vCommon[i]);



                // This new combination is this new graph + all possible combinations of subsequent graphs retreived by
                // the recursive call
                kNewComb.push_back( kTotal );

                kRes.push_back( kNewComb );

            }



            iC++;
            // Generate all combinations, selecting iK-1 from the vertex set
            // -1 is because we've already selected one common vertex
        //} while( next_combination(kComb.begin(), kComb.begin() + (iK-1), kComb.end()) );



    }

    return make_pair(kRes, true);
}



double GraphHandler::VertexRank(Vertex V) const
{
    //map< Vertex, double > rank;
    return out_degree(V, *m_RankG);
    //Graph& G = *m_RankG;

    // Loop through all the vertices in the full graph
    // For each, update their rank with the given equation until some threshold is found

    /*
    // Calculate page rank
    double dPRank = 1.0;
    double c = 0.1;

    typedef graph_traits<Graph>::out_edge_iterator Edg_It;
    Edg_It first, last;

    // Add all the edges to the new vertex
    // ------------------------------------
    for (tie(first, last) = out_edges(V, G); first != last; ++first)
    {
        int iTarget = get(vertex_name, G, target(*first, G));

    }*/
}

/*
  Based on algorithm RaRe from the paper "FINDING COMMUNITIES BY CLUSTERING A GRAPH
INTO OVERLAPPING SUBGRAPHS"

Consider a connected component C at iteration i. If C has more vertices than our
maximum desired core size max, we remove a set Ri of vertices, where | Ri |= t


  */
void GraphHandler::ClusterComponent(Graph H)
{
    int iMaxClusterSize = 3;
    int iMinClusterSize = 2;


    int t=1;
    if( num_vertices(H) > iMaxClusterSize )
    {
        // {vi} are the t highest ranking nodes in H
        // 1. Compute the rankings of the nodes in H

        Iter_t vi, vi_end;
        tie(vi, vi_end) = vertices(H);

        vector< pair<unsigned int, double> > v;
        v.reserve(num_vertices(H));
        m_RankG = &H;
        for( tie(vi, vi_end) = vertices(H); vi != vi_end; ++vi )
            v.push_back(pair<unsigned int, double>(get(vertex_name, H, *vi), VertexRank(*vi)));



        sort( v.begin(), v.end(), sort_pred<Graph>() );

        //---------------------------
        set<unsigned int>& R = m_RSet;
        // 2. R <- R U {vi}; H <- H \{vi}
        for(int i=0; i < t; i++)
        {
            //R.insert( v[i].first );
            R.insert( v[i].first );

            //H <- H \{vi}
            cout << "Removing " << m_IDTypes[v[i].first] << endl;
            Vertex vTmp = FindVertex(v[i].first, H).second;
            clear_vertex(vTmp, H); // Remove all edges from the vertex
            remove_vertex(vTmp, H); // Remove the vertex from the graph
        }


        static int iDebugC = 0;
        stringstream sFName;
        sFName << "debug/debuggraph" << iDebugC << ".dot";
        SaveGraph(sFName.str(), H);
        iDebugC++;
        cout << "Saving " << sFName.str() << endl;

        // {Fi} are connected components in H
        // for all {Fi}, ClusterComponent(Fi)
        std::vector<int> component(num_vertices(H));
        int num = connected_components(H, &component[0]);

        //for (int i = 0; i != component.size(); ++i)
        //      cout << "Vertex " << i <<" is in component " << component[i] << endl;



        for(int i=0; i < num; i++)
        {
            vector<Vertex> Fi;
            for(int j=0; j != component.size(); ++j)
                if( component[j] == i )
                    Fi.push_back(j);

            // Fi is a list of vertices for the i:th component
            Graph k = InducedSubgraph(Fi, H);
            cout << "Calling cluster component" << endl;
            ClusterComponent(k);
        }


    }
    else// if( iMinClusterSize <= num_vertices(H) && num_vertices(H) <= iMaxClusterSize )
    {
        // Mark H as cluster core
        m_ClusterCores.push_back(H);
    }
}

GraphHandler::Graph GraphHandler::InducedSubgraph(vector<Vertex> V, const Graph& G)
{
    Graph SG;

    map<Vertex, Vertex> kMap;

    // Add the vertex set
    for(int i=0; i < V.size(); i++)
    {
        Vertex vNew = add_vertex(SG);
        put(vertex_name, SG, vNew,get(vertex_name, G, V[i]));
        kMap[V[i]] = vNew;
    }

    // Add edges
    graph_traits<Graph>::edge_iterator firstE, lastE;
    for (tie(firstE, lastE) = edges(G); firstE != lastE; ++firstE)
    {
        bool bSFound = kMap.find(source(*firstE, G)) != kMap.end();
        bool bTFound = kMap.find(target(*firstE, G)) != kMap.end();
        if( bSFound && bTFound )
        {
            // Add the edge, both are in induced subgraph
            add_edge(kMap[source(*firstE, G)], kMap[target(*firstE, G)], SG);
        }
    }

    return SG;
}
/*
  Based on algorithm RaRe from the paper "FINDING COMMUNITIES BY CLUSTERING A GRAPH
INTO OVERLAPPING SUBGRAPHS"
  */
vector<GraphHandler::Graph> GraphHandler::RareAlgorithm(Graph& G)
{
    vector<Graph> C; // Clusters

    m_RSet.clear();
    ClusterComponent(G);



    // The initial clusters {Ci} are the cluster cores
    for( int i=0; i < m_ClusterCores.size(); i++ )
        C.push_back(m_ClusterCores[i]);

    // for all v in R
    for( set<unsigned int>::iterator it = m_RSet.begin(); it != m_RSet.end(); ++it )
    {
        // For all clusters Ci do:
        for( int i=0; i < C.size(); i++ )
        {
            // Add v to Ci if v is adjacent to Ci or W(v U Ci)>W(Ci)

            // 1. Check if v is adjacent to Ci
            graph_traits<Graph>::adjacency_iterator v1, v2;
            Vertex vR = FindVertex(*it, G).second;
            bool bFound = false;
            for( tie(v1, v2) = adjacent_vertices(vR, G); (v1 != v2) && !bFound; ++v1 )
                // Check if we can find any vertex in Ci from vR, if so they are adjacent
                if( FindVertex(get(vertex_name, G, *v1) , C[i] ).first )
                    bFound = true;

            bFound = false;
            if( bFound )
            {
                // Add v to Ci
                Vertex vNew = add_vertex(C[i]);
                put(vertex_name, C[i], vNew, get(vertex_name, G, vR));

                cout << "Adding vertex: " << m_IDTypes[get(vertex_name, G, vR)] << endl;

                // Add all the edges as well

                typedef graph_traits<Graph>::out_edge_iterator Edg_It;
                Edg_It first, last;


                // Add all the edges to the new vertex
                // ------------------------------------
                for (tie(first, last) = out_edges(vR, G); first != last; ++first)
                {
                    int iTarget = get(vertex_name, G, target(*first, G));
                    Vertex vTarget;
                    bool bFound = false;
                    tie(bFound, vTarget) = FindVertex(iTarget, C[i]);
                    // If the target vertex exist (out of R) in C[i] and the edge does not exist
                    if( bFound )
                    {
                        cout << "Found target" << endl;
                        pair<Edge, bool> kRet = edge( vNew, vTarget, C[i]);

                        // Add if it doesn't exist
                        if( !kRet.second )
                        {
                            cout << "Adding edge" << endl;
                            kRet = add_edge(vNew, vTarget, C[i]);
                        }
                    }
                }
                // ------------------------------------
            }
        }
    }

    return C;
}

void GraphHandler::SaveMetisGraph(string sFilename, const Graph& G)
{
    ofstream out(sFilename.c_str(), ios_base::out);

    out << num_vertices(G) << " " << num_edges(G) << endl;

    Iter_t vi, vi_end;
    // For each vertex
    for (tie(vi, vi_end) = vertices(G); vi != vi_end; ++vi)
    {
        //out << "v " << j << " " << get(vertex_name, G, *vi) << endl;

        // Print all adjacent vertices
        graph_traits<Graph>::adjacency_iterator v1, v2, v3;
        for( tie(v1, v2) = adjacent_vertices(*vi, G); v1 != v2; ++v1 )
        {

                v3 = v1; v3++;
                if( v3 != v2 ) // If there is another adjacent vertex after this one
                    //out << get(vertex_name, G, *v1) << " ";
                    out << (*v1 +1) << " ";
                else
                    out << (*v1 +1); // Last one
                    //out << get(vertex_name, G, *v1); // Last one

        }
        out << endl;
    }

    out.close();
}

vector<GraphHandler::Graph> GraphHandler::LoadMetisPartition(string sFilename, const Graph& G)
{
    vector<Graph> kPartition;
    vector< vector<Vertex> > vList;
    map< Vertex, int> kPartMap;

    int iNumPart = -1; // Number of partitions

    ifstream in(sFilename.c_str());
    Iter_t vi, vi_end;
    for (tie(vi, vi_end) = vertices(G); vi != vi_end; ++vi)
    {
        string sTmp;
        getline(in, sTmp, '\n');

        int iPart = lexical_cast<int>(sTmp);
        //vList
        kPartMap[*vi] = iPart;

        if( iPart > iNumPart )
            iNumPart = iPart;
    }

    iNumPart++;
    cout << "iNumPart: " << iNumPart << endl;
    vList.resize(iNumPart);
    // Now create the graphs
    for(int i=0; i < iNumPart; i++)
    {
        for (tie(vi, vi_end) = vertices(G); vi != vi_end; ++vi)
        {
            if( kPartMap[*vi] == i )
            {
                vList[i].push_back(*vi);
            }
        }
        kPartition.push_back( InducedSubgraph(vList[i], G) );
    }

    in.close();
    return kPartition;

}

vector< pair<Cluster*, double> > GraphHandler::GetPatterns(const Graph& G)
{
    set< pair<Cluster*, double> > pPatterns;
    vector<unsigned int> iIDVector;
    map<Cluster*, unsigned int> kClusterMap;

    int iMinC = 2;

    // Get all closed frequent subgraphs
    for(int i=0; i < m_pAllclusters.size(); i++)
    {
        if( IsSubgraphIsomorph(m_pAllclusters[i]->kMeta, G)
                //&& num_vertices(m_pAllclusters[i]->kMeta) >= iMinC )
                && num_vertices(m_pAllclusters[i]->kMeta) >= iMinC )
        {
            // So this specific subgraph is contained in the graph.
            // Check: 1. If this specific subgraph is contained in any other pattern already added, don't add
            // If we pass check 1, we add it and remove all of those patterns which are contained in this specific subgraph
            // Pair consist of (Cluster*, frequency)
            //pPatterns.push_back( make_pair(m_pAllclusters[i], m_pAllclusters[i]->pGraphs.size()) );

            // Check 1.
            bool b1 = true;
            for( set< pair<Cluster*, double> >::iterator it = pPatterns.begin(); it != pPatterns.end() && b1; ++it )
            {
                if( IsSubgraphIsomorph(m_pAllclusters[i]->kMeta, it->first->kMeta) )
                    b1 = false;
            }

            if( b1 )
            {
                // Add it and remove all of those patterns which are contained in this specific subgraph
                for( set< pair<Cluster*, double> >::iterator it = pPatterns.begin(); it != pPatterns.end(); ++it )
                {
                    // If the pattern *it is contained in the one we're adding, remove it (since the one we're adding is bigger)
                    if( IsSubgraphIsomorph(it->first->kMeta, m_pAllclusters[i]->kMeta) )
                    {
                        pPatterns.erase(it);
                    }
                }

                pPatterns.insert( make_pair(m_pAllclusters[i], m_pAllclusters[i]->pGraphs.size()) );
                kClusterMap[m_pAllclusters[i]] = i;
            }
        }
    }

    //cout << "Mined closed graphs: " << pPatterns.size() << endl;

    vector< pair<Cluster*, double> > kRes;
    kRes.reserve(pPatterns.size());
    for( set< pair<Cluster*, double> >::iterator it = pPatterns.begin(); it != pPatterns.end(); ++it )
    {
        kRes.push_back(*it);
        iIDVector.push_back(kClusterMap[it->first]);
    }

    m_iIDVector = iIDVector;


    return kRes;
}

vector<GraphHandler::Graph> GraphHandler::DecomposeGraph2(const Graph& G)
{
    vector<Graph> kRes;
    //vector< pair<Cluster*, double> > pPatterns;
    int iMinC = min(2, (int)num_vertices(G)); // The minimum size we're considering

    set< pair<Cluster*, double> > pPatterns;

    // Get all closed frequent subgraphs
    for(int i=0; i < m_pAllclusters.size(); i++)
    {
        if( IsSubgraphIsomorph(m_pAllclusters[i]->kMeta, G)
                //&& num_vertices(m_pAllclusters[i]->kMeta) >= iMinC )
                && num_vertices(m_pAllclusters[i]->kMeta) >= iMinC )
        {
            // So this specific subgraph is contained in the graph.
            // Check: 1. If this specific subgraph is contained in any other pattern already added, don't add
            // If we pass check 1, we add it and remove all of those patterns which are contained in this specific subgraph
            // Pair consist of (Cluster*, frequency)
            //pPatterns.push_back( make_pair(m_pAllclusters[i], m_pAllclusters[i]->pGraphs.size()) );

            // Check 1.
            bool b1 = true;
            for( set< pair<Cluster*, double> >::iterator it = pPatterns.begin(); it != pPatterns.end() && b1; ++it )
            {
                if( IsSubgraphIsomorph(m_pAllclusters[i]->kMeta, it->first->kMeta) )
                    b1 = false;
            }

            if( b1 )
            {
                // Add it and remove all of those patterns which are contained in this specific subgraph
                for( set< pair<Cluster*, double> >::iterator it = pPatterns.begin(); it != pPatterns.end(); ++it )
                {
                    // If the pattern *it is contained in the one we're adding, remove it (since the one we're adding is bigger)
                    if( IsSubgraphIsomorph(it->first->kMeta, m_pAllclusters[i]->kMeta) )
                        pPatterns.erase(it);
                }

                pPatterns.insert( make_pair(m_pAllclusters[i], m_pAllclusters[i]->pGraphs.size()) );
            }
        }
    }

    //cout << "Mined " << pPatterns.size() << " closed frequent subgraphs." << endl;

    kRes.reserve(pPatterns.size());
    for( set< pair<Cluster*, double> >::iterator it = pPatterns.begin(); it != pPatterns.end(); ++it )
    {
        kRes.push_back(it->first->kMeta);
    }

    if( pPatterns.size() == 0 )
        kRes.push_back(G);

    //vector<Graph> GiTmp;
    Graph kHyp = GraphUnion(kRes);
    SaveGraph("graphunion.dot", kHyp);
    Graph kDiff = GetDifferenceGraph(G, kHyp);
    SaveGraph("diffg.dot", kDiff);

    // Now get the components of kDiff
    std::vector<int> component(num_vertices(kDiff));
    int num = connected_components(kDiff, &component[0]);

    // For each of these components, we will add them to GiTmp
    // First we need to form the induced subgraphs of them, since they are just
    // assignement of each vertex to a component number
    for(int i=0; i < num; i++)
    {
        vector<Vertex> Fi;
        for(int j=0; j != component.size(); ++j)
            if( component[j] == i )
                Fi.push_back(j);

        // Fi is a list of vertices for the i:th component
        kRes.push_back(InducedSubgraph(Fi, kDiff));
    }

    //GetDifferenceGraph
    //kFinal.push_back(GiTmp);

    // Make a matrix of the pairwise graph edit distances between the patterns
    // Select the highest ones

    return kRes;
}

vector<GraphHandler::Graph> GraphHandler::DecomposeGraph3(const Graph& G)
{
    vector< pair<Cluster*, double> > pPatterns;
    pPatterns = GetPatterns(G);


    vector<Graph> C;

    if( pPatterns.size() == 0 )
    {
        C.push_back(G);
        return C;
    }

    // 2. Choose the frequent subgraph of the highest frequency

    sort( pPatterns.begin(), pPatterns.end(), sort_pred<Cluster*>() );

    //---------------------
    for(int i=0; i < pPatterns.size(); i++)
    {
        C.push_back( pPatterns[i].first->kMeta );
    }
    return C;

    //---------------------

    // Since it is sorted according to ascending order, the highest is at the end.
    C.push_back( pPatterns.back().first->kMeta );
    pPatterns.pop_back();
    //--------------------------------------------------

    Graph kUnion;
    while(num_vertices(kUnion) < num_vertices(G))
    {
        kUnion = GraphUnion(C);
        bool bFound;
        bFound = false;

        for(int i=0; i < C.size(); i++)
        {
            pair<bool, Graph*> Gres;

            // Check if we can find another frequent subgraph to connect to G_i
            // The criterion is that they share atleast one vertex and that adding it would make the union graph
            // "bigger"
            if( (Gres = FindCommonFreqSubgraph(pPatterns, C[i], kUnion, C)).first )
            {
                C.push_back(*Gres.second);
                bFound = true;
                break;
            }
        }
        if(!bFound)
        {

            //-----------------------------------------
            // Did not find any
            //Graph kHyp = GraphUnion(Gi);

            if( num_vertices(G) > num_vertices(kUnion) )
            {

            Graph kDiff = GetDifferenceGraph(G, kUnion);
            //SaveGraph("diffg.dot", kDiff);
            // Now get the components of kDiff
            std::vector<int> component(num_vertices(kDiff));
            int num = connected_components(kDiff, &component[0]);

            // For each of these components, we will add them to GiTmp
            // First we need to form the induced subgraphs of them, since they are just
            // assignement of each vertex to a component number
            for(int i=0; i < num; i++)
            {
                vector<Vertex> Fi;
                for(int j=0; j != component.size(); ++j)
                    if( component[j] == i )
                        Fi.push_back(j);

                // Fi is a list of vertices for the i:th component
                C.push_back(InducedSubgraph(Fi, kDiff));
            }

            }
            //-----------------------------------------
            return C;
        }
    }

    return C;

}

vector<GraphHandler::Graph> GraphHandler::DecomposeGraph(const Graph& G)
{
    // 1. Check which frequent subgraphs are occuring in G
    // --------------------------------------------
    vector<Graph> kRes;
    vector< pair<Cluster*, double> > pPatterns;
    int iMinC = min(2, (int)num_vertices(G)); // The minimum size we're considering


    // Atleast 200
    pPatterns.reserve(200);

    //cout << "Num vertices G:" << num_vertices(G) << endl;

/*
    for(int i=0; i < m_pAllclusters.size(); i++)
    {
        if( IsSubgraphIsomorph(m_pAllclusters[i]->kMeta, G)
                //&& num_vertices(m_pAllclusters[i]->kMeta) >= iMinC )
                && num_vertices(m_pAllclusters[i]->kMeta) >= iMinC )
            // Pair consist of (Cluster*, frequency)
            pPatterns.push_back( make_pair(m_pAllclusters[i], m_pAllclusters[i]->pGraphs.size()) );
    }*/
    //cout << "Matched: " << pPatterns.size() << endl;

    pPatterns = GetPatterns(G);




    if( pPatterns.size() == 0 )
    {
        //GetOptimalParentOnline
        //GetOptimalParentOnline(G);
        kRes.push_back(G);
        return kRes;
    }

    // Select on random
    //--------------------------- RANDOM
    /*
    vector<Graph> kResTmp;

    int iRSize = 1 + min(rand() % (int)pPatterns.size(),10);
    if( iRSize == pPatterns.size() )
        iRSize = pPatterns.size() - 1;
    random_shuffle(pPatterns.begin(), pPatterns.end());
    //kResTmp.resize(iRSize);
    //vector< pair<Cluster*, double> > pPatterns;
    for(int i=0; i < iRSize; i++)
        kResTmp.push_back(pPatterns[i].first->kMeta);
    return kResTmp;*/
    //---------------------------

    // --------------------------------------------
    // 2. Choose the frequent subgraph of the highest frequency

    sort( pPatterns.begin(), pPatterns.end(), sort_pred<Cluster*>() );

    // Since it is sorted according to ascending order, the highest is at the end.
    kRes.push_back( pPatterns.back().first->kMeta );
    pPatterns.pop_back();

    // Now find the most frequent subgraph which has some vertex in common with kRes (if it exists)
    //pair<bool, Graph*> G2 = FindCommonFreqSubgraph(pPatterns, kRes[0], kRes[0]);

    //if( G2.first )
    //    kRes.push_back(*G2.second);
   // else
  //      return kRes;



    vector<bool> bFinished;
    list< vector<Graph> > kList;
    list< vector<Graph> > kFinished;

    // Each of the solutions will start with kRes, but then we may connect to 1...k different
    // of the already added kRes

    //kRes.clear();
    vector< vector<GraphHandler::Graph> > kFinal = DecomposeGraphInternal(kRes, G, pPatterns);
    for(int i=0; i < kFinal.size(); i++)
        kFinal[i].push_back(kRes.back());

    /*
    while(true)
    {
        // For each of the different total solutions, we may form

        for( list< vector<Graph> >::iterator it=kList.begin(); it!=kList.end(); ++it )
        {
            //Graph kUnion = GraphUnion(kRes);
            Graph kUnion = GraphUnion(*it);

            // Mark one solution as done
            if( num_vertices(kUnion) >= num_vertices(G) )
            //if( IsSubgraphIsomorph(kUnion, G) )
            {
                kFinished.push_back(*it);
                kList.erase(it);
                //break;
            }
        }

        bool bFound = false;

        // Used to store the graphs temp
        vector<Graph> kTmp;
        kTmp.clear();
        // 3. Now, given a list of graphs, we might connect a new subgraph G_i+1
        // to any of the already existing subgraphs
        for(int i=0; i < kRes.size(); i++)
        {
            //vector< vector<GraphHandler::Graph> > kTmp = DecomposeGraphInternal(kRes, kRes[i], pPatterns);
            // We then combine the results
            pair<bool, Graph*> Gi;

            if( (Gi = FindCommonFreqSubgraph(pPatterns, kRes[i], kUnion)).first )
            {
                // They do have some vertex in common, which means its a valid addition
                //DecomposeGraphInternal()
                //kRes.push_back(*Gi.second);
                kTmp.push_back(*Gi.second);
                bFound = true;
                break;
            }
        }

        // If we cannot find any more
        if( !bFound )
            break;

        for(int i=0; i < kTmp.size(); i++)
            kRes.push_back(kTmp[i]);

    }*/

    // Find best
    vector< Graph >* pBest = 0;
    int iBestDist = 9999;
    int iChosenIndex = -1;

    // TODO: Make it choose the most highly correlated "chain"
    for(int i=0; i < kFinal.size(); i++)
    {
        int iDist;


        Graph kUnion = GraphUnion(kFinal[i]);
        //SaveGraph("graphunion.dot", kUnion);
        /*
        if( (iDist = EditDistance(kUnion, G)) < iBestDist )
        {
            pBest = &kFinal[i];
            iBestDist = iDist;
        }*/
        if( (iDist = kFinal[i].size()) < iBestDist && IsSubgraphIsomorph(G, kUnion) )
        {
            iChosenIndex = i;
            pBest = &kFinal[i];
            iBestDist = iDist;
        }
    }

    //cout << "Total kFinal.size(): " << kFinal.size() << "; Chosen: " << iChosenIndex << endl;
    pBest = 0;
    if( pBest == 0 )
    {
        iBestDist = 9999;
        pBest = 0;
        for(int i=0; i < kFinal.size(); i++)
        {
            int iDist;
            Graph kUnion = GraphUnion(kFinal[i]);
            //SaveGraph("graphunion.dot", kUnion);
            if( (iDist = EditDistance(kUnion, G)) < iBestDist && kFinal[i].size() != 1 )
            {
                pBest = &kFinal[i];
                iBestDist = iDist;
            }
        }
    }

    if( kFinal.size() != 0 )
    {
        if( pBest != 0 )
        return *pBest;//kFinal[0];
        else
            return kFinal[0];
    }
    else
        return kRes;
}



vector< vector<GraphHandler::Graph> > GraphHandler::DecomposeGraphInternal(vector<Graph> Gi,const Graph& G,
           vector< pair<Cluster*, double> > pPatterns)
{
    // We have a vector of subgraphs thus far

    Graph kUnion = GraphUnion(Gi);
    if( num_vertices(kUnion) >= num_vertices(G) )
        return vector< vector<GraphHandler::Graph> >();

    vector< vector<GraphHandler::Graph> > kFinal;

    // 3. Given a list of graphs, we might connect a new subgraph G_i+1
    // to any of the already existing subgraphs
    //for(int i=0; i < Gi.size(); i++)
    for( vector<Graph>::iterator it = Gi.begin(); it != Gi.end(); ++it )
    {
        //kFinal.push_back(Gi);
        pair<bool, Graph*> Gres;

        // Check if we can find another frequent subgraph to connect to G_i
        // The criterion is that they share atleast one vertex and that adding it would make the union graph
        // "bigger"
        if( (Gres = FindCommonFreqSubgraph(pPatterns, *it, kUnion, Gi)).first )
        {
            vector<Graph> GiTmp;
            //vector<Graph> GiTmp = Gi;
            GiTmp.push_back(*Gres.second);
            //SaveGraph("tmpgraph.dot", *Gres.second);
            // They do have some vertex in common, which means its a valid addition

            vector<Graph> GiTmp3 = Gi;
            GiTmp3.push_back(*Gres.second);
            // Keep generating all subsequent possibilities as well
            //vector< vector<GraphHandler::Graph> > kRes = DecomposeGraphInternal(GiTmp, G, pPatterns);
            vector< vector<GraphHandler::Graph> > kRes = DecomposeGraphInternal(GiTmp3, G, pPatterns);

            if( kRes.size() != 0 )
            {
                // Add kRes to the final list combined with kFinal[i]
                //for(int j=0; j < kRes.size(); j++)
                for(int j=0; j < 1; j++)
                {
                    if( kRes[j].size() != 0 )
                    {
                    vector<Graph> GiTmp2 = GiTmp;
                    //GiTmp
                    GiTmp2.resize(GiTmp2.size() + kRes[j].size());
                    copy(kRes[j].begin(), kRes[j].end(), GiTmp2.begin() + GiTmp.size());
                    kFinal.push_back(GiTmp2);
                    return kFinal;
                    }
                }
            }
            else
            {
                // TODO: Check if G=GraphUnion(Gi), if not we need to add some clusters which might not be FS.

                //GetDifferenceGraph
                kFinal.push_back(GiTmp);
                return kFinal;
            }

        }
        else
        {

        }
    }

/*
    vector<Graph> GiTmp;

    Graph kHyp = GraphUnion(Gi);
    //SaveGraph("graphunion.dot", kHyp);

    if( num_vertices(G) > num_vertices(kHyp) )
    {

    Graph kDiff = GetDifferenceGraph(G, kHyp);
    //SaveGraph("diffg.dot", kDiff);
    // Now get the components of kDiff
    std::vector<int> component(num_vertices(kDiff));
    int num = connected_components(kDiff, &component[0]);

    // For each of these components, we will add them to GiTmp
    // First we need to form the induced subgraphs of them, since they are just
    // assignement of each vertex to a component number
    for(int i=0; i < num; i++)
    {
        vector<Vertex> Fi;
        for(int j=0; j != component.size(); ++j)
            if( component[j] == i )
                Fi.push_back(j);

        // Fi is a list of vertices for the i:th component
        GiTmp.push_back(InducedSubgraph(Fi, kDiff));
    }

    }

    kFinal.push_back(GiTmp);*/

    return kFinal;

}


pair<bool, GraphHandler::Graph*> GraphHandler::FindCommonFreqSubgraph(vector< pair<Cluster*, double> >& pPatterns, const Graph& G, const Graph& kUnion, vector<Graph>& Gi) const
{
    Graph* kG;
    bool bFound = false;

    //double dBest = -10.0;
    //Graph* pBest = 0;

    for(int i=pPatterns.size()-1; i >= 0; --i )
    {
        // Check so that we have atleast some vertex in common
        if( !HasVertexCommon(G, pPatterns[i].first->kMeta, 0) )
        {
            // Check so that its not already included
            bool bValid = true;
            for(int j=0; j < Gi.size() && bValid; j++)
            {
                if( EditDistance(pPatterns[i].first->kMeta, Gi[j]) == 0 )
                    bValid = false;
                // If this pattern which is added is already contained in some other added Gi, don't add it
                if( IsSubgraphIsomorph(pPatterns[i].first->kMeta, Gi[j]) )
                    bValid = false;
            }



            if( !IsSubgraphIsomorph(pPatterns[i].first->kMeta, kUnion) && bValid )
            {
                kG = &pPatterns[i].first->kMeta;
                bFound = true;
                // Remove all those patterns which are contained in this one
                pPatterns.erase(pPatterns.begin() + i);


                for(int j=0; j < pPatterns.size(); j++)
                    if( IsSubgraphIsomorph(pPatterns[j].first->kMeta, *kG) )
                        pPatterns.erase(pPatterns.begin() + j);


/*
                for(int j=0; j < Gi.size(); j++)
                {
                    if( IsSubgraphIsomorph(Gi[j], *kG) )
                        Gi.erase(Gi.begin() + j);
                }*/


                break;
            }
        }
    }

    return make_pair<bool, Graph*>(bFound, kG);
}

/*
Creates a file and saves all those vertices which should be blacklisted according to the threshold
  */
void GraphHandler::SaveBlacklist(string sFilename, int iThreshold)
{
    set<string> kBlacklist = InternalGenerateBlacklist(iThreshold);

    ofstream out(sFilename.c_str(), ios_base::out);

    for( set<string>::iterator it = kBlacklist.begin(); it != kBlacklist.end(); ++it )
        out << *it << endl;

    out.close();
}

/*
Loads a file with names of the vertices to blacklist and then performs the blacklist on all graphs
  */
void GraphHandler::LoadBlacklist(string sFilename)
{
    set<string> kBlacklist;

    ifstream in(sFilename.c_str(), ios_base::out);

    //for( set<string>::iterator it = kBlacklist.begin(); it != kBlacklist.end(); ++it )
    //    out << *it << endl;

    while( !in.eof() )
    {
        string sTmp;
        getline(in, sTmp, '\n');
        kBlacklist.insert(sTmp);
    }

    in.close();

    cout << "Starting blacklist..." << endl;
    for( set<string>::iterator it = kBlacklist.begin(); it != kBlacklist.end(); ++it )
    {
        set<string> tmpSet;
        tmpSet.insert(*it);
        cout << "Blacklisting: " << *it << endl;
        Blacklist(tmpSet);
        if( m_Types.find(*it) != m_Types.end() )
            m_Types.erase(m_Types.find(*it));
    }

    cout << "Blacklisted " << kBlacklist.size() << " vertices. Remaining: " << m_Types.size() << endl;

}

/*
    // Go through all the different category types, then check in the whole cluster database to see
     if the category type is "frequent enough"
    // We might for example have the criterion that the category type is existing in atleast one FS(frequent subgraph) of size >= 2.

  */
set<string> GraphHandler::InternalGenerateBlacklist(int iThreshold) const
{
    set<string> kBlacklist;


    for( set<string>::const_iterator it = m_Types.begin(); it != m_Types.end(); ++it )
    {
        bool bFreq = false;
        Graph gTmp;
        Vertex vTmp = add_vertex(gTmp);
        map<string, int>::const_iterator iTmpIt = m_TypesID.find(*it);
        //put(vertex_name, gTmp, vTmp, m_TypesID[*it]);
        put(vertex_name, gTmp, vTmp, (*iTmpIt).second);

        cout << "Checking: " << *it << endl;

        // So, now check if this type is "frequent enough"
        int iFreq = 0;
        for(int i=0; i < m_pAllclusters.size() && !bFreq; i++)
            if( HasVertexCommon(m_pAllclusters[i]->kMeta, gTmp, 1) && m_pAllclusters[i]->pGraphs.size() >= 2 )
                iFreq++;

        bFreq = (iFreq > iThreshold);

        if( !bFreq )
        {
            kBlacklist.insert(*it);
            cout << *it << endl;
        }

    }

    return kBlacklist;
}



void GraphHandler::BlacklistOnCriterion(int iThreshold)
{
    set<string> kBlacklist = InternalGenerateBlacklist(iThreshold);

    cout << "Starting blacklist..." << endl;
    for( set<string>::iterator it = kBlacklist.begin(); it != kBlacklist.end(); ++it )
    {
        set<string> tmpSet;
        tmpSet.insert(*it);
        cout << "Blacklisting: " << *it << endl;
        Blacklist(tmpSet);
        if( m_Types.find(*it) != m_Types.end() )
            m_Types.erase(m_Types.find(*it));
    }

    cout << "Blacklisted " << kBlacklist.size() << " vertices. Remaining: " << m_Types.size() << endl;

}

GraphHandler::Graph GraphHandler::FormResultOfVertexEditOp(const Graph& GP, pair<Vertex, Vertex> vEditOp)
{
    Graph kRes = GP;
    Vertex vExisting = FindVertex(vEditOp.first, kRes).second;
    Vertex vNew = add_vertex(kRes);
    put(vertex_name, kRes, vNew, vEditOp.second);



    // Now add the edge as well
    add_edge(vExisting, vNew, kRes);

    //-------------------
    kEditG = Graph();
    Vertex v1 = add_vertex(kEditG);
    Vertex v2 = add_vertex(kEditG);
    put(vertex_name, kEditG, v1, vEditOp.second);
    put(vertex_name, kEditG, v2, get(vertex_name, GP, vExisting));

    add_edge(v1, v2, kEditG);
    //-------------------
    return kRes;
}

GraphHandler::Graph GraphHandler::FormResultOfEdgeEditOp(const Graph& GP, pair<unsigned int, unsigned int> eEditOp)
{
    Graph kRes = GP;
    Vertex v1 = FindVertex(eEditOp.first, kRes).second;
    Vertex v2 = FindVertex(eEditOp.second, kRes).second;

    // Now add the edge as well
    add_edge(v1, v2, kRes);


    //-------------------
    kEditG = Graph();
    Vertex vu = add_vertex(kEditG);
    Vertex vv = add_vertex(kEditG);

    put(vertex_name, kEditG, vu, get(vertex_name, kRes, v1));
    put(vertex_name, kEditG, vv, get(vertex_name, kRes, v2));

    add_edge(vu, vv, kEditG);
    //-------------------

    return kRes;
}

pair<GraphHandler::Graph, int> GraphHandler::GetOptimalParentOnline(const Graph& GP,
  const vector<Graph>& kVec, int iIndex, const Graph& G, const Graph* pUnion)
{
    // Lets go through all the graphs in the graph database, and see what the optimal
    // parent might be.

    // The "optimal" next thing would be the one edit distance away operation
    // such that it is most frequent among all those instance in which we find the graph G

    // Given that we've only seen G as the partial graph, what action would we chose?

    // 1. Given G as the partial graph in some G_i in the graph data base
    // Enumerate all possible next actions

    // Save all edge addition actions in one set
    // Save all vertex addition actions in another

    // Sum everything up


    int iTotalGraphsExamined = 0;

    typedef pair<Vertex, Vertex> VEditOp;
    map< VEditOp, int > kSum;

    typedef pair<unsigned int, unsigned int> EEditOp;
    map< EEditOp, int > kESum;


    for(int i=0; i < m_Graphs.size(); i++)
    {
        if( IsSubgraphIsomorph(GP, m_Graphs[i]) && i != iExcludeGraph )
        {
            iTotalGraphsExamined++;
            GraphExploration GE(m_Graphs[i], GP);
            GraphExploration::EditOpIterator it(GE);
            // Go through all possible vertex addition operations
            while( !it.isDoneVertex() )
            {
                pair<Vertex, Vertex> vEdit = it.NextVertexOp();
                if( kSum.find(vEdit) == kSum.end() )
                    kSum[vEdit] = 1;
                else
                    kSum[vEdit]++;
            }

            // Go through all possible edge addition operations
            while( !it.isDoneEdge() )
            {
                EEditOp eEdit = it.NextEdgeOp();
                if( kESum.find(eEdit) == kESum.end() )
                    kESum[eEdit] = 1;
                else
                    kESum[eEdit]++;
            }
        }
    }

    Graph kTotalUnion;

    if( pUnion == 0 )
        kTotalUnion = GraphUnion(kVec);
    else
        kTotalUnion = *pUnion;

    //cout << "Examined " << iTotalGraphsExamined << " different graphs" << endl;
    VEditOp kBest;
    EEditOp kEBest;
    int iVBest = -1; // The value of the best vertex op
    int iEBest = -1; // The value of the best edge op

    //ofstream out("theta.dat", ios_base::out);

    //out << (int)kSum.size() << " " << (int)kESum.size() << endl;

    int iTotalEditSum = 0;
    for( map< VEditOp, int >::iterator it=kSum.begin(); it != kSum.end(); ++it )
        iTotalEditSum += it->second;

    for( map< EEditOp, int >::iterator it=kESum.begin(); it != kESum.end(); ++it )
        iTotalEditSum += it->second;

    vector< pair<string, double> > editStats;

    //-------------------------------------------------------
    // Check which was the best vertex edit op, and save how good it was
    for( map< VEditOp, int >::iterator it=kSum.begin(); it != kSum.end(); ++it )
    {
        Graph vOpRes = FormResultOfVertexEditOp(GP, it->first);

        //out <<  m_IDTypes[it->first.first] << ", " << m_IDTypes[it->first.second] << '\t' << ((double)it->second / (double)iTotalEditSum)  << endl;
        stringstream sTmp; sTmp << m_IDTypes[it->first.first] << ", " << m_IDTypes[it->first.second];
        editStats.push_back( make_pair(sTmp.str(), (double)it->second / (double)iTotalEditSum));

        // vEditOp.second, ID of the new vertex to connect to
        // vEditOp.first, ID of the existing vertex
        if( it->second > iVBest && !IsSubgraphIsomorph(vOpRes, kTotalUnion) )
        {
            //---------------------//
            // Make sure that the new graph union contains some new information
            vector<Graph> kTmpVec;
            if( pUnion == 0 )
            {
            for(int j=0; j < kVec.size(); j++)
                if( j != iIndex )
                    kTmpVec.push_back(kVec[j]);
            }
            else
                kTmpVec.push_back(kTotalUnion);
            // Also add what would be the result of this addition
            kTmpVec.push_back(vOpRes);
            Graph kUnion = GraphUnion(kTmpVec);
            if( !GraphEqual(G, kUnion) )
            {
                iVBest = it->second;
                kBest = it->first;
            }
            //---------------------//
        }
    }
    //-------------------------------------------------------



    //-------------------------------------------------------
    // Check which was the best edge edit op, and save how good it was
    for( map< EEditOp, int >::iterator it=kESum.begin(); it != kESum.end(); ++it )
    {
        // This check is to ensure that adding this edge actually makes a difference in the "complete" graph
        Graph vOpRes = FormResultOfEdgeEditOp(GP, it->first);

        //out <<  m_IDTypes[it->first.first] << ", " << m_IDTypes[it->first.second] << '\t' << ((double)it->second / (double)iTotalEditSum)  << endl;
        stringstream sTmp; sTmp << m_IDTypes[it->first.first] << ", " << m_IDTypes[it->first.second];
        editStats.push_back( make_pair(sTmp.str(), (double)it->second / (double)iTotalEditSum));
        if( it->second > iEBest && !IsSubgraphIsomorph(vOpRes, kTotalUnion) )
        {
            //---------------------//
            // Make sure that the new graph union contains some new information
            vector<Graph> kTmpVec;
            for(int j=0; j < kVec.size(); j++)
                if( j != iIndex )
                    kTmpVec.push_back(kVec[j]);
            // Also add what would be the result of this addition
            kTmpVec.push_back(vOpRes); Graph kUnion = GraphUnion(kTmpVec);
            if( !GraphEqual(G, kUnion) )
            {
                // It does contribute with "new information", record it
                iEBest = it->second;
                kEBest = it->first;
            }
            //---------------------//
        }
    }
    //-------------------------------------------------------

    sort(editStats.begin(), editStats.end(), sort_pred< string >() );
    //for(int i=0; i < editStats.size(); i++)
    //    out <<  editStats[i].first << '\t' << editStats[i].second  << endl;

    //out.close();

    if( iVBest >= iEBest )
    {

        //cout << "Optimal parent online, u: " << m_IDTypes[kBest.first] << ", v: " << m_IDTypes[kBest.second] << endl;
        //cout << "Support: " << iVBest << endl;
        return make_pair(FormResultOfVertexEditOp(GP, kBest), iVBest);
    }
    else
    {
        // TODO: Fix edge edit ops as well
        //cout << "Optimal parent online, edge op; u: " << m_IDTypes[kEBest.first] << ", v: " << m_IDTypes[kEBest.second] << endl;
        //cout << "Support: " << iEBest << endl;
        return make_pair(FormResultOfEdgeEditOp(GP, kEBest), iEBest);
    }



}

int GraphHandler::GetTotalNumVertices() const
{
    int iV = 0;
    BOOST_FOREACH( const Graph& G, m_Graphs )
        iV += num_vertices(G);
    return iV;
}

int GraphHandler::GetTotalNumEdges() const
{
    int iV = 0;
    BOOST_FOREACH( const Graph& G, m_Graphs )
        iV += num_edges(G);
    return iV;
}

int GraphHandler::GetTotalNumLabels() const
{
    return m_Types.size();
}

int GraphHandler::GetTotalNumGraphs() const
{
    return m_Graphs.size();
}

void GraphHandler::AssertUniqueLabels()
{
    int iG = 0;
    BOOST_FOREACH( Graph& G, m_Graphs )
    {
        std::pair<Iter_t, Iter_t> vp;
        set<unsigned int> added;
        for (vp = vertices(G); vp.first != vp.second; ++vp.first)
        {
            int iID = get(vertex_name, G, *vp.first);
            assert( added.find(iID) == added.end() );
        }
        iG++;
    }
}

void GraphHandler::RemoveDisconnectedGraphs()
{
    vector<Graph> kNewGraphs;
    BOOST_FOREACH( Graph& G, m_Graphs )
    {
        std::vector<int> component(num_vertices(G));
        int num = connected_components(G, &component[0]);
        if( num == 1 )
            kNewGraphs.push_back(G);
    }
    cout << "Removed " << m_Graphs.size() - kNewGraphs.size() << " disconnected graphs" << endl;
    m_Graphs = kNewGraphs;
}

GraphHandler::Graph GraphHandler::BruteForcePredict(const Graph* pPartialGraph)
{

    m_bCorrectPrediction = false;
    if( pPartialGraph == 0 )
        pPartialGraph = &m_kPartialG;
    Graph G = *pPartialGraph;

    int iBest = -1;
    Graph kBest;

    // For each vertex, check what is the most optimal vertex to add
    std::pair<Iter_t, Iter_t> vp;
    vector<Graph> vTmp;
    int iIndex = -1;
    for (vp = vertices(G); vp.first != vp.second; ++vp.first)
    {
        Graph gTmp;
        Vertex vNew = add_vertex(gTmp);
        put(vertex_name, gTmp, vNew, get(vertex_name, G, *vp.first));
        vTmp.push_back(gTmp);
    }

    // For each vertex, check what is the most optimal vertex to add
    for(int i=0; i < vTmp.size(); i++)
    {
        //-----------
        pair<Graph, int> kOTmp = GetOptimalParentOnline(vTmp[i], vTmp, i, G, pPartialGraph);
        //-----------

        if( kOTmp.second > iBest )
        {
            iIndex = i;
            //kHyp.push_back(*(pParent.first));
            iBest = kOTmp.second;
            kBest = kOTmp.first;
        }
    }

    //SaveGraph("kbest.dot", kBest);
    //SaveGraph("tmp.dot", vTmp[iIndex]);

    if( IsSubgraphIsomorph(kBest, *m_pEGraph) && !IsSubgraphIsomorph(kBest, *pPartialGraph) && iBest != -1 )
        m_bCorrectPrediction = true;

    vector<Graph> kTmpVec;
    /*
    for(int j=0; j < vTmp.size(); j++)
        if( j != iIndex )
            kTmpVec.push_back(vTmp[j]);*/
    kTmpVec.push_back(*pPartialGraph);
    // Also add what would be the result of this addition
    kTmpVec.push_back(kBest);
    Graph kUnion = GraphUnion(kTmpVec);



    return kUnion;
}

// Returns the number of times G has been observed in the graph database
int GraphHandler::GetGraphSupportOnline(const Graph& G) const
{
    int iRes = 0;
    for(int i=0; i < m_Graphs.size(); i++)
        if( IsSubgraphIsomorph(G, m_Graphs[i]) && i != iExcludeGraph )
            iRes++;
    return iRes;
}

double GraphHandler::GetGraphCorrelation(int iG1Supp, int iG2Supp, int iG1G2Supp) const
{
    if( iG1Supp == 0 || iG1Supp == 1 || iG2Supp == 0 || iG2Supp == 1 )
        return 0.0;
    double dSupG1G2 = (double)iG1G2Supp / (double)m_Graphs.size();
    double dSupG1 = (double)iG1Supp / (double)m_Graphs.size();
    double dSupG2 = (double)iG2Supp / (double)m_Graphs.size();
        return (dSupG1G2 - dSupG1*dSupG2) / sqrt(dSupG1*dSupG2*(1-dSupG1)*(1-dSupG2));
}

double GraphHandler::GetGraphCorrelation(const Graph& G1, const Graph& G2) const
{
    int iG1 = GetGraphSupportOnline(G1);
    int iG2 = GetGraphSupportOnline(G2);
    int iG1G2Supp = JointGraphSupp(G1, G2);
    //cout << "iG1: " << iG1 << ", " << "iG2: " << iG2 << ", G1G2Supp: " << iG1G2Supp << endl;
    return GetGraphCorrelation(iG1, iG2, iG1G2Supp);
}

int GraphHandler::JointGraphSupp(const Graph& G1, const Graph& G2) const
{
    //static vector<Graph*> kTmp;
    int iRes = 0;
    for(int i=0; i < m_Graphs.size(); i++)
    {
        bool b1 = false, b2 = false;
        if( IsSubgraphIsomorph(G1, m_Graphs[i]) )
            b1 = true;
        if( IsSubgraphIsomorph(G2, m_Graphs[i]) )
            b2 = true;
        if( b1 && b2 )
            iRes++;
    }

    return iRes;
}

/*
double GraphHandler::GetGraphCorrelation(const Graph& G1, const Graph& G2)
{

}
*/

pair<GraphHandler::Graph, double> GraphHandler::GetOptimalParentOnline2(const Graph& GP,
  const vector<Graph>& kVec, int iIndex, const Graph& G, const Graph* pUnion)
{
    // Lets go through all the graphs in the graph database, and see what the optimal
    // parent might be.

    // The "optimal" next thing would be the one edit distance away operation
    // such that it is most frequent among all those instance in which we find the graph G

    // Given that we've only seen G as the partial graph, what action would we chose?

    // 1. Given G as the partial graph in some G_i in the graph data base
    // Enumerate all possible next actions

    // Save all edge addition actions in one set
    // Save all vertex addition actions in another

    // Sum everything up


    int iTotalGraphsExamined = 0;

    typedef pair<Vertex, Vertex> VEditOp;
    map< VEditOp, int > kSum;

    typedef pair<unsigned int, unsigned int> EEditOp;
    map< EEditOp, int > kESum;


    /*
    double dSupG1 = (double)(m_pAllclusters[i])->pGraphs.size();
    double dSupG2 = (double)(m_pAllclusters[j])->pGraphs.size();

    double dSupG1G2 = JointGraphSupp(m_pAllclusters[i], m_pAllclusters[j]);// / (double)m_Graphs.size();

    return (dSupG1G2 - dSupG1*dSupG2) / sqrt(dSupG1*dSupG2*(1-dSupG1)*(1-dSupG2));
      */

    int iSuppG1 = 0;
    // kSum[vEdit] gives iSuppG1G2 where G2 is the resulting graph of vEdit
    // What is G2? The number of times G2 has been observed in total


    // (dSupG1G2 - dSupG1*dSupG2) / sqrt(dSupG1*dSupG2*(1-dSupG1)*(1-dSupG2));

    for(int i=0; i < m_Graphs.size(); i++)
    {
        if( IsSubgraphIsomorph(GP, m_Graphs[i]) && i != iExcludeGraph )
        {
            iTotalGraphsExamined++; iSuppG1++;
            GraphExploration GE(m_Graphs[i], GP);
            GraphExploration::EditOpIterator it(GE);
            // Go through all possible vertex addition operations
            while( !it.isDoneVertex() )
            {
                pair<Vertex, Vertex> vEdit = it.NextVertexOp();
                if( kSum.find(vEdit) == kSum.end() )
                    kSum[vEdit] = 1;
                else
                    kSum[vEdit]++;
            }

            // Go through all possible edge addition operations
            while( !it.isDoneEdge() )
            {
                EEditOp eEdit = it.NextEdgeOp();
                if( kESum.find(eEdit) == kESum.end() )
                    kESum[eEdit] = 1;
                else
                    kESum[eEdit]++;
            }
        }
    }

    Graph kTotalUnion;

    if( pUnion == 0 )
        kTotalUnion = GraphUnion(kVec);
    else
        kTotalUnion = *pUnion;

    //cout << "Examined " << iTotalGraphsExamined << " different graphs" << endl;
    VEditOp kBest;
    EEditOp kEBest;
    int iVBest = -1; // The value of the best vertex op
    int iEBest = -1; // The value of the best edge op
    double dVBest = -2.0;
    double dEBest = -2.0;

    //-------------------------------------------------------
    // Check which was the best vertex edit op, and save how good it was
    for( map< VEditOp, int >::iterator it=kSum.begin(); it != kSum.end(); ++it )
    {
        Graph vOpRes = FormResultOfVertexEditOp(GP, it->first);
        int iG2Supp = GetGraphSupportOnline(vOpRes);

        // it->second now has the joint support
        //GetGraphCorrelation(iSuppG1, iG2Supp, it->second );
        double dCorr = GetGraphCorrelation(iSuppG1, iG2Supp, it->second );
        //if( it->second > iVBest && !IsSubgraphIsomorph(vOpRes, kTotalUnion) )
        if( dCorr > dVBest && !IsSubgraphIsomorph(vOpRes, kTotalUnion) )
        {
            //---------------------//
            // Make sure that the new graph union contains some new information
            vector<Graph> kTmpVec;
            if( pUnion == 0 )
            {
            for(int j=0; j < kVec.size(); j++)
                if( j != iIndex )
                    kTmpVec.push_back(kVec[j]);
            }
            else
                kTmpVec.push_back(kTotalUnion);
            // Also add what would be the result of this addition
            kTmpVec.push_back(vOpRes);
            Graph kUnion = GraphUnion(kTmpVec);
            if( !GraphEqual(G, kUnion) )
            {
                dVBest = dCorr;
                kBest = it->first;
            }
            //---------------------//
        }
    }
    //-------------------------------------------------------

    //-------------------------------------------------------
    // Check which was the best edge edit op, and save how good it was
    for( map< EEditOp, int >::iterator it=kESum.begin(); it != kESum.end(); ++it )
    {
        // This check is to ensure that adding this edge actually makes a difference in the "complete" graph
        Graph vOpRes = FormResultOfEdgeEditOp(GP, it->first);
        int iG2Supp = GetGraphSupportOnline(vOpRes);
        double dCorr = GetGraphCorrelation(iSuppG1, iG2Supp, it->second );
        //if( it->second > iEBest && !IsSubgraphIsomorph(vOpRes, kTotalUnion) )
        if( dCorr > dEBest && !IsSubgraphIsomorph(vOpRes, kTotalUnion) )
        {
            //---------------------//
            // Make sure that the new graph union contains some new information
            vector<Graph> kTmpVec;
            for(int j=0; j < kVec.size(); j++)
                if( j != iIndex )
                    kTmpVec.push_back(kVec[j]);
            // Also add what would be the result of this addition
            kTmpVec.push_back(vOpRes); Graph kUnion = GraphUnion(kTmpVec);
            if( !GraphEqual(G, kUnion) )
            {
                // It does contribute with "new information", record it
                //iEBest = it->second;
                dEBest = dCorr;
                kEBest = it->first;
            }
            //---------------------//
        }
    }
    //-------------------------------------------------------

    //if( iVBest >= iEBest )
    if( dVBest >= dEBest )
    {

        //cout << "Optimal parent online, u: " << m_IDTypes[kBest.first] << ", v: " << m_IDTypes[kBest.second] << endl;
        //cout << "Support: " << iVBest << endl;
        return make_pair(FormResultOfVertexEditOp(GP, kBest), dVBest);
    }
    else
    {
        // TODO: Fix edge edit ops as well
        //cout << "Optimal parent online, edge op; u: " << m_IDTypes[kEBest.first] << ", v: " << m_IDTypes[kEBest.second] << endl;
        //cout << "Support: " << iEBest << endl;
        return make_pair(FormResultOfEdgeEditOp(GP, kEBest), dEBest);
    }



}

void GraphHandler::GenerateFrequentEditMaps(const vector<GSpanStruct>& vFreqSubgraphs)
{
    vector<EditOperationMap>& kEditMaps = m_kEditMaps;
    for(unsigned int i=0; i < vFreqSubgraphs.size(); i++)
    {
        cout << i << "/" << vFreqSubgraphs.size() << endl;
        kEditMaps.push_back(EditOperationMap(vFreqSubgraphs[i].kGraphIDs, vFreqSubgraphs[i].G));
    }

}
